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