r语言进行go富集分析_生信实操|如何利用R语言进行GSEA分析

看图说话栏目曾介绍过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分析，先要对