r语言进行go富集分析

看图说话栏目曾介绍过GSEA的原理(看图说话|GSEA分析--教你解锁高级的富集分析)，今天我们来看一下如何利用R语言进行GSEA分析。

如果你有RNA-seq的数据，就可以这样做，先把数据整成这样，一共两列，一列是SYMBOL，一列是foldChange。

然后输入代码：

> setwd("E:\\ 20201214GSEA分析怎么做")
> #if (!requireNamespace("BiocManager", quietly = TRUE))
> #  install.packages("BiocManager")
> #BiocManager::install("clusterProfiler")
> library(clusterProfiler)#主要用到clusterProfiler包
> df = read.table("easy_input.txt",header = T,as.is = T)
> head(df)#查看前面几行
    SYMBOL  foldChange
1      A2M -0.50361737
2     NAT1 -3.25012047
3     NAT2 -0.28492113
4 SERPINA3 -0.63137249
5    AADAC -0.18896102
6     AAMP -0.03605823
> dim(df)#数据总共几行几列
[1] 12444     2
>

可以看出，输入的数据总共12444个基因，当然，这并不仅仅是是差异基因。
在分析之前，还是要把SYMBOL转换成ENTREZ ID，之前讲过，ENTREZ ID在染色体上是有编号的，一个编号对应一个基因，错不了，但是如果是gene name，很容易出错，另外也不好计算。
这里要用到clusterProfiler包里面的bitr函数，这个函数很方便的就可以转换ENTREZID和SYMBOL。

> df.id#转换的列是df数据框中的SYMBOL列
+             fromType = "SYMBOL",#需要转换ID类型
+             toType = "ENTREZID",#转换成的ID类型
+             OrgDb = "org.Hs.eg.db")#对应的物种，小鼠的是org.Mm.eg.db
'select()' returned 1:many mapping between keys and columns
Warning message:
In bitr(df$SYMBOL, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db") :
  1.28% of input gene IDs are fail to map
> head(df.id)
    SYMBOL ENTREZID
1      A2M        2
2     NAT1        9
3     NAT2       10
4  SERPINA3       12
5    AADAC       13
6     AAMP       14
> dim(df.id)
[1] 12287     2
>

代码说有1.28%没有map，所以，刚才是12444个基因，现在就变成了了12287个基因。
接下来，将最开始的文件和转成ENTREZID的文件合并，用到merge函数。

> #让基因名、ENTREZID、foldchange对应起来
> easy.df"SYMBOL",all=F)
> head(easy.df)
  SYMBOL  foldChange ENTREZID
1   A1CF -0.13606162    29974
2    A2M -0.50361737        2
3 A4GALT  0.22796773    53947
4  A4GNT -0.08822687    51146
5   AAAS -0.12161476     8086
6   AACS  0.23240224    65985

为了做GSEA分析，先要对所有的基因进行排序，就是大家常看到的那个蝴蝶图，这里用foldchange作为量化的标准，将基因从大到小的排序。

> #按照foldchange排序
> sortdf$foldChange, decreasing = T),]#降序排序
> head(sortdf)
     SYMBOL foldChange ENTREZID
6398   MMP1   4.572613     4312
1724  CDC45   4.514594     8318
7042    NMU   4.418218    10874
1731  CDCA8   4.144075    55143
6190  MCM10   3.876258    55388
1706  CDC20   3.677857      991
> gene.expr = sortdf$foldChange#把foldchange按照从大到小提取出来
> head(gene.expr)
[1] 4.572613 4.514594 4.418218 4.144075 3.876258 3.677857
> names(gene.expr) $ENTREZID#给上面提取的foldchange对应上ENTREZID
> head(gene.expr)
    4312     8318    10874    55143    55388      991 
4.572613 4.514594 4.418218 4.144075 3.876258 3.677857

有了这个ENTREZID和foldchange信息的gene.expr之后，就可以使用clusterProfiler进行GSEA分析了，进行GSEA分析。
这里就很简单，GSEA就一句代码：

kk "hsa")

是不是出奇的简单。

> kk "hsa")
preparing geneSet collections...
GSEA analysis...
leading edge analysis...
done
Warning messages:
1: In preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam,  :
  There are ties in the preranked stats (0.02% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.
2: In serialize(data, node$con) :
  'package:stats' may not be available when loading
3: In serialize(data, node$con) :
  'package:stats' may not be available when loading
4: In serialize(data, node$con) :
  'package:stats' may not be available when loading
5: In serialize(data, node$con) :
  'package:stats' may not be available when loading
6: In fgseaMultilevel(...) :
  For some pathways, in reality P-values are less than 1e-10. You can set the `eps` argument to zero for better estimation.
> head(kk)
               ID                  Description setSize enrichmentScore      NES
hsa04110 hsa04110                   Cell cycle     114       0.6698529 2.740968
hsa03050 hsa03050                   Proteasome      43       0.7099227 2.425653
hsa05169 hsa05169 Epstein-Barr virus infection     193       0.4337991 1.932303
hsa04657 hsa04657      IL-17 signaling pathway      85       0.5624088 2.172820
hsa03030 hsa03030              DNA replication      33       0.7233029 2.309149
hsa01230 hsa01230  Biosynthesis of amino acids      62       0.5891945 2.178783
               pvalue     p.adjust      qvalues rank                   leading_edge
hsa04110 1.000000e-10 3.240000e-08 2.431579e-08 1210 tags=40%, list=10%, signal=37%
hsa03050 1.358588e-08 2.200913e-06 1.651757e-06 2469 tags=65%, list=20%, signal=52%
hsa05169 8.736724e-08 9.435662e-06 7.081345e-06 2770 tags=39%, list=23%, signal=31%
hsa04657 1.413224e-07 1.144711e-05 8.590914e-06 2830 tags=49%, list=23%, signal=38%
hsa03030 2.526900e-07 1.637431e-05 1.228871e-05 1867 tags=64%, list=15%, signal=54%
hsa01230 1.176431e-06 6.352728e-05 4.767642e-05 2868 tags=55%, list=23%, signal=42%
                                                                                                                                                                                                                                                                                                                                                                               core_enrichment
hsa04110                                                                                                                                                         8318/991/9133/890/983/4085/7272/1111/891/4174/9232/4171/993/990/5347/701/9700/898/23594/4998/9134/4175/4173/10926/6502/994/699/4609/5111/1869/1029/8317/4176/2810/3066/1871/1031/9088/995/1019/4172/5885/11200/7027/1875/7534
hsa03050                                                                                                                                                                                                                                       5688/5709/5698/5693/3458/5713/11047/5721/5691/5685/5690/5684/5686/5695/10213/23198/7979/5699/5714/5702/5708/5692/5704/5683/5694/5718/51371/5682
hsa05169 3627/890/6890/9636/898/9134/6502/6772/3126/3112/4609/917/5709/1869/3654/919/915/4067/4938/864/4940/5713/5336/11047/3066/54205/1871/578/1019/637/916/3383/4939/10213/23586/4793/5603/7979/7128/6891/930/5714/3452/6850/5702/4794/7124/3569/7097/5708/2208/8772/3119/5704/7186/5971/3135/1380/958/5610/4792/10018/8819/3134/10379/9641/1147/5718/6300/3109/811/5606/2923/3108/5707/1432
hsa04657                                                                                                                                                                 4312/6280/6279/6278/3627/2921/6364/8061/4318/3576/3934/6347/727897/1051/6354/3458/6361/6374/2919/9618/5603/7128/1994/7124/3569/8772/5743/7186/3596/6356/5594/4792/9641/1147/2932/6300/5597/27190/1432/7184/64806/3326
hsa03030                                                                                                                                                                                                                                                                           4174/4171/4175/4173/2237/5984/5111/10535/1763/5427/23649/4176/5982/5557/5558/4172/5424/5983/5425/54107/6119
hsa01230                                                                                                                                                                                                             29968/26227/875/445/5214/440/65263/6472/7086/5723/3418/7167/586/2597/5230/2023/5223/5831/6888/50/5315/3419/10993/2805/22934/3421/5634/3417/5232/221823/2027/5211/384/5832

返回的kk文件很复杂,总共48行，11列。

> dim(kk)
[1] 48 11

1，代表的KEGG中的信号通路，比如hsa04110就是KEGG中的Cell cycle的信号通路
2，对信号通路的描述，也就是信号通路的名字
3，该信号通路的基因个数
4，很熟悉的富集分数，就是ES
5，这是标准化以后的ES，全称是normalized enrichment score
6，P值
7，矫正以后的p值
8，Q值，也有的写FDR q-val(false discovery rate)，表示错误发现率
9，这里的rank就是ES在顶点时候，那个最幸运基因的位置，或者说是排名
10，这里的leading_edge有点复杂，对富集贡献最大的基因成员,即领头亚集，用于定义Leading-edge subset的参数有：Tags，List，Signal，对于一个基因集而言，定义其中对Enrichment score贡献最大的基因为核心基因，也称之为leading edge subset，也就是ES顶点的之前的基因。tags表示核心基因占该基因集基因总数的比例，而list表示核心基因占所有基因总数的比例，signal利用这两个指标计算得到，公式如下：

N代表所有基因的数目，Nh代表该基因集下的基因总数。对于一个基因集而言，当核心基因的数目和该基因集下的基因总数相同，signal取值最大，当该基因集的基因数目和所有基因数目接近时，signal的取值接近于0。当然，我们希望的是signal越大越好。
11，这个core_enrichment就是主要富集的基因，下面的数字其实是ENTREZID，就是一个一个基因，换句话说，就是做GSEA分析的基因列表中，Hit到该目的通路的基因列表。

接下里，需要按照ES的大小，给通路排序，然后保存起来，后面用来可视化。

#按照enrichment score从高到低排序，便于查看富集通路
sortkkT),]
head(sortkk)
dim(sortkk)
write.table(sortkk,"gsea_output2.txt",sep = "\t",quote = F,col.names = T,row.names = F)

通过查看gsea_output2.txt,选择你想画的pathway，记下通路的term ID
开始画图,不需要另外准备输入文件，clusterProfiler直接可以用gseaplot出图：

gseaplot(kk, "hsa04510")#很简单的代码是吧

可以改变参数，更改线条的颜色。

gseaplot(kk, "hsa04510", color.line='steelblue')

在这里，我们看到ES是负值，ES是负值，就是表示，功能集中在尾部，就是表示实验组低表达，对照组高表达。
这里用到clusterProfiler画图，如果觉得不够美观，还可以用enrichplot的包来画图

#另一种出图方式是模拟Broad institute的GSEA软件的图，可以用gseaplot2出图：
library(enrichplot)
gseaplot2(kk, "hsa04510")
gseaplot2(kk, "hsa04510", color = "firebrick", rel_heights=c(1, .2, .6))#改变更多参数，为了美观

是不是感觉有点Nature的味道了。
它还有更外一个功能，就是支持同时展示多个pathways的结果。

gseaplot2(kk, row.names(sortkk)[1:4])

上图用的是ES排名前4个画图，也可以用你自己感兴趣的通路画图，自己在刚才保存的txt文件里挑选就行。

paths "hsa04510", "hsa04512", "hsa04974", "hsa05410")
gseaplot2(kk, paths)

这里的GSEA分析其实由三个图构成，GSEA分析的runningES折线图，还有下面基因的位置图，还有所谓的蝴蝶图。如果不想同时展示，还可以通过subplots改变。

gseaplot2(kk, paths, subplots=1)#只要第一个图
gseaplot2(kk, paths, subplots=1:2)#只要第一和第二个图
gseaplot2(kk, paths, subplots=c(1,3))#只要第一和第三个图

如果想把p值标上去，也是可以的。

gseaplot2(kk, paths, color = colorspace::rainbow_hcl(4), pvalue_table = TRUE)

到目前为止，基本上用KEGG库做注释的GSEA分析就做的差不多了，还想了解更多的，可能需要对clusterProfiler ，enrichplot两个包好好研读一下，说不定有惊喜。
最后甩出一个代码：

gseaplot2(kk,#数据
          row.names(sortkk)[39],#画那一列的信号通路
          title = "Prion disease",#标题
          base_size = 15,#字体大小
          color = "green",#线条的颜色
          pvalue_table = TRUE,#加不加p值
          ES_geom="line")#是用线，还是用d点

注：此推文未经许可禁止转载！ 阅读推荐： 文献详解|又一篇7+的纯生信分析文章，带你解析肿瘤预后模型研究套路工具篇 | TransCirc：一个预测分析circRNA翻译的数据库

期刊推荐|最快4周接收，这本不收版面费的期刊今年将突破5分，值得投稿

工具篇|批量下载文献法宝①——关键词篇