clusterProfiler包进行KEGG,GO,GSEA富集分析

本地的KEGG分析参考文章：KEGG数据库使用及通路分析教程,GO参考文章：FunRich数据库：一个主要用于基因和蛋白质的功能富集以及相互作用网络分析的独立的软件工具，当然该工具不止可以进行富集分析，具体去看文章吧。

这里就先介绍一下本地GSEA分析

我们前文说过，GSEA分析是基于表达矩阵的。所以我们首先得有一个基因表达矩阵。除此以外，我们至少还需要一个表型文件，其实就是表达数据的分组信息而已。

对于表达数据的文件格式有4种：

• GCT：基因簇文本文件格式（* .gct）
• RES：ExpRESsion（带P和A调用）文件格式（* .res）
• PCL：斯坦福cDNA文件格式（* .pcl）
• TXT：表达式数据集的文本文件格式（* .txt）

可能我们最常用的就是txt格式的数据文件。我们这里以txt的文件讲解。格式如下：

第一列是基因名称，也就是symbol，第二列是描述性信息，可以都是na，其他列都是不同样本标准化后的表达数据。

对于另外3个格式文件，大家可以参考其他公众号的文章，这里给大家分享一篇：手把手教你GSEA分析。也不用你自己再去搜索了，我再写一遍也没有意义。

第二个文件是表型文件，格式是cls格式文件，此文件可以根据样品对应的表型按照下面文件格式自己制作，但是，文件后缀必须为*. cls，且用tab或者空格分割的文本文件。

第一行是3个数字，第一个数值数我们样本数，第二个是我们分组个数，我们这里只有Normal 和Tumor组所以是2，第三个是1，这是固定不变的，写上就行。

第二行是以#号开头，后面就是我们的分组名称。

第三行就是我们437个样本是Normal还是Tumor，顺序和表达数据文件中的一样。

我们以之前上传的TCGA数据库33个Project的RNA-Seq转录组数据为例，选择TCGA-COAD进行分析，TCGA转录组数据处理方式，参考文章：TCGA数据库：RNA-Seq数据的下载与处理。

rm(list = ls())
options(stringsAsFactors = F)
load("F:/TCGA/HTSeq-FPKM/Rdata/data/TCGA-COAD-Exp.Rdata")
protExp <- transomeData[["proteinCodingExpData"]][["Exp"]]
protExp <- data.matrix(protExp) %>% as.data.frame()
protExp <- protExp[apply(t(protExp),2,var)!=0,]
GSEAExp <- data.frame(NAME = rownames(protExp),Description = "na",protExp,check.names = F)
GSEAExp[1:5,1:4]
write.table(GSEAExp,file= "GSEAExp.txt",sep = "\t",row.names = F,quote = F)
NormalNum <- transomeData[["proteinCodingExpData"]][["sampleNum"]][["NormalNum"]]
TumorNum <-transomeData[["proteinCodingExpData"]][["sampleNum"]][["TumorNum"]]
cat(file = "GSEAExp.cls",paste((ncol(GSEAExp)-2),2,1),"\n#Normal Tumor\n",rep("Normal",NormalNum),rep("Tumor",TumorNum))

上面代码是处理之前准备的数据，每个人的数据可能不一样，所以代码会不同，总之，我们处理后会得到上面2个作为GSEA分析的输入文件。

一.本地软件GSEA分析

本地分析的话，需要下载软件。或者谷歌/百度搜一下。

官网：https://www.gsea-msigdb.org/gsea/index.jsp

如果第一次使用需要点击“Click here” 注册，注册过后登录，只需要在文本框内输入注册邮箱，点击login登录即可。

根据你的系统选择下载安装：

安装后打开软件：

load data我们的数据。GSEAExp.txt是我们的表示数据，GSEAExp.cls是分组信息。

加载数据后，如果我们的数据没有问题，就会提示成功：

在Run GSEA窗口设置参数。

Expression dataset（表达文件）：选择上一步上传的GSEAExp.txt文件

Gene sets database （功能基因集数据库）：GSEA包含了MSigDB数据库中的功能基因集，可以从中选择感兴趣的通路、癌症标记、转录因子数据库等。我们在前面文章:为什么选择GSEA分析？和KEGG和GO分析有什么区别？中就介绍了这些数据集，当然，这个数据集我们可以自己准备，多数情况下，我们是选择数据库给我们定义好的数据集，所以直接用就好了。我们这里选择：c2.cp.kegg.v6.2.symbols.gmt

Number of permutations（扰动／随机次数）：通常设置1000，此参数不可过小。默认就行。

Phenotypes labels(样品表型分类文件)：选择上一步上传的表型cls文件

我们选择一个就行，2个分组的话，其实没多大区别。

Collapse dataset to gene symbols：我们选择No_Collapse

Permutation type(扰动类型)：通常选择phenotype，如果样品数目较少可以选择gene_set。

Chip platform(芯片类型)：如果表达gct文件的第一列为芯片探针id则此处需要选择对应的芯片平台，如果是基因symbol则无需选择。

然后点击Run,左下角就可以看见程序开始运行了,时间可能很久。

运行结束后显示成功，点击我们就可以进入一个页面。就是我们的分析结果。

文件也在本地找的到，就是在我们电脑的用户名文件下的gsea_home/output中，比如：

file:///C:/Users/Administrator/gsea_home/output/jun11/my_analysis.Gsea.1591844474693/index.html

部分结果如下：

经过各项参数筛选后剩下176个基因集的条目，其中有124个条目在Tumor组中上调（ES值为正）；有35个基因集的False Discovery Rate （FDR）小于0.25，是有意义的。在未校正的p<0.01下，有14个基因集富集显著，在未校正的p<0.05下，有21个基因集富集显著；富集结果浏览，可点击Snapshot展示了ES绝对值最大的20个基因集的图。点击enrichment results in html format查看详细的网页形式的富集分析结果。点击 enrichment results in excel查看详细的Excel表格形式的富集结果。点击Guide to的网页链接查看官方对结果解读的指导文档。

对于Enrichment in phenotype: Normal (39 samples)解释和上面一样。