春节后的第一次写学习笔记，还是要提醒大家没事别聚集、别聚集、别聚集！重要的事情说三遍！
前几天看到公众号里写了GSEA的文章，虽然我之前跟着很多网上的教程走过一遍GSEA分析，但是对于GSEA的结果图并不知道怎么看，所以也想学习一下这方面的知识。
在跟着教程走代码之前，想先了解一下什么是GSEA，以及这个东西是干嘛用的。
（一）什么是GSEA?
GSEA: Gene Set Enrichment Analysis
翻译应该不难吧： 基因集富集分析
（二）上面说的集是什么意思？
参考文章：GSEA入门------原理
在以前的实验中发表的数据或表达谱上共表达的基因信息数据集合（瞧这专业的说法，我看了好几遍才看明白）。说的通俗一点：就是某一个通路（相关的所有基因的总和）。
（三）GSEA分析的是什么东东？
我们在利用DESeq2，或者edgeR，或者limma进行差异分析之后，会得到一个列表，里面有你的所有样品的差异基因。我们对这个列表里的基因表达水平进行一个排序，从高到低。GSEA分析就是确定一个基因集S（或者一个通路）的成员是否倾向于出现在列表L的顶部(或底部)，如果聚集在顶部/底部，该基因集与表型分类区分相关。GSEA是一种无阈值方法，可根据其差异表达等级或其他分数对所有基因进行分析，无需事先进行基因过滤（学习基因通路富集分析软件GSEA）。
（四）GSEA原理是什么？
实在是不想看一堆公式。。。所以就把教程里（刘小泽学习GSEA）的一段话再说的明白点：
假设：某一个通路的全部基因在排序后的差异基因列表中随机分布。
但是：如果我们看到它们“意外地”出现在基因列表的某一端（从图上看是在某一侧形成一个峰），那么就可以计算显著性来看看富集程度如何。如果富集结果显著，那么就拒绝原假设，认为这个通路的基因在我们的基因列表中富集，并且看到富集分数（是个发paper的契机啊）。
（五）上面提到的富集分数是什么意思？
富集分数：Enrichment Score，简称ES。
从我们排序的差异基因列表里，从高到低一个一个来看：遇到一个基因，属于我们要的基因集S时，ES就会加分；如果遇到的这个基因不属于我们要的基因集，ES就会减分。来一张图，更直观一些：
我们通常看到的GSEA的图是长这样的：
（六）实际操作：方法一：R语言
这里我用的是DESeq2得到的差异基因列表，这个列表是怎么得来的，请参考：TCGA的差异基因分析。如果你有自己的差异基因列表更好。这个差异基因的列表可以是DESeq2，也可以是edgeR或者limma分析得到的。
#GSEA analysis steps
#得到差异分析结果
> geneList <- DESeq2_DEG$log2FoldChange
#调用需要的R包
> library(stringr)
> library(clusterProfiler)
> library(org.Hs.eg.db)
先看一下我的这个DESeq2差异基因列表，可以看到基因名称都是Ensemble ID：
我们需要把基因名称先转换成Entrez ID:
#如果是Ensemble ID，并且如果还带着版本号，需要先去除版本号
> names(geneList) <- str_split(rownames(DESeq2_DEG) ,
pattern = '\\.',simplify = T)[,1]
#下面用clusterProfiler包进行基因名称的转换：
> geneList_tr <- bitr(names(geneList),
fromType = "ENSEMBL",
toType = c("ENTREZID","SYMBOL"),
OrgDb = org.Hs.eg.db)
> new_list <- data.frame(ENSEMBL=names(geneList),
logFC = as.numeric(geneList))
> new_list <- merge(new_list, geneList_tr, by = "ENSEMBL")
> geneList <- new_list$logFC
> names(geneList) <- geneList_tr$ENTREZID
# 最后从大到小排序，得到一个字符串
> geneList <- sort(geneList,decreasing = T)
#进行GO分析
> go_result <- gseGO(geneList = geneList,
ont = "BP",
OrgDb = org.Hs.eg.db,
minGSSize = 10,
maxGSSize = 500,
pvalueCutoff = 1)
#将其中的基因名变成symbol ID
> go_result <- setReadable(go_result, OrgDb = org.Hs.eg.db)
#还可以直接点击查看，只需要转成数据框
> go_result_df <- as.data.frame(go_result)
#如果对结果的某个通路感兴趣，可以作图
> gseaplot(go_result, geneSetID = c("GO:0010638"))
#进行KEGG分析
> kk <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 10,
maxGSSize = 500,
pvalueCutoff = 1,
verbose = FALSE)
> go_result_kk <- as.data.frame(kk)
> gseaplot(kk, geneSetID = "hsa05202")
（七）实际操作：GSEA软件(这个方法麻烦得很啊，还是R语言方便快捷，但做出来的图漂亮，有闲工夫的可以试试）
（1）软件下载
GSEA的官网：https://www.gsea-msigdb.org/gsea/index.jsp
这里以windows系统为例。安装软件就不说了，so easy。
运行一下吧，界面长这样：
（2）准备需要的文件
先别着急用，为什么我说这个软件很烦躁呢，因为它不是啥格式的文件都识别的。
首先，是你的基因列表。这里的基因列表请注意！不是差异基因列表，是你所有的基因！！！（参考文章：Question: GSEA error ）可以先打开看一下我们的表（这个表是怎么来的，可以参照利用R包TCGAbiolinks进行各种数据下载去下载RNA-seq的FPKM结果来进行练习）：
然而这样的表是不能够被GSEA软件识别的。根据官网的说明：Dataformats:
The TXT format is a tab delimited file format that describes an expression dataset. It is organized as follows:
The first line contains the labels Name and Description followed by the identifiers for each sample in the dataset. NOTE: The Description column is intended to be optional, but there is currently a bug such that it is treated as required. We hope to fix this in a future release. If you have no descriptions available, a value of NA will suffice.
注意注意注意！！！这里的示范列表里比我们多出了一列DESCRIPTION。这一列是必须有的。你要自己添加进去。如果没有什么可描述的，就写上na。如果这里你不注意，之后软件会一直报错，报到你怀疑人生！！！
你的基因列表还需要符合下面的格式：
第一行: Name(tab)Description(tab)(sample 1 name)(tab)(sample 2 name) (tab) ... (sample N name)
NOTE：必须是制表符分隔！别的符号统统不可以！
那么怎么让它变成制表符分隔呢？请看：
在excel里点击“导出”，再点击“更改文件类型”：
选择“文本文件（制表符分隔）”，最后点击“另存为”就好了。要确保是.txt后缀的文件格式。
接下来还需要一个文件，就是你的分组信息。在前面的文章里，提到了如何把TCGA数据库里的样品分组（TCGA的差异基因分析）。那么这个分组信息的文件格式也是有讲究的，先看一下官网上给出的示范：
那么根据官网的例子，我也整理了一下我的分组信息：
NOTE:这里需要注意的是，这个文件你需要保存成.cls后缀的文件。
（3）万事俱备，开始分析
有了必须的文件，就可以进行分析了。
首先，上传你的文件（method1选择你的文件）：
如果你的文件格式没问题，这里会提示你成功上传了2个文件。
然后点击左边一栏的"Run GSEA":
然后就是一系列的选择参数了：
然后选择Gene sets database，点右边的...会出来一个选择框：
是不是觉得好多？选择困难症要犯了。。。我好南。。。和我一样的童鞋可以看这个官网：https://www.gsea-msigdb.org/gsea/msigdb，这里对每一个database都进行了说明：
但看完感觉还是不知道选哪个怎么办？我问了一下实验室的Steven，他说比如说你的样品是control和knockout，你想看看这两者的区别，那么就选择H，因为这里面是一些经典的生物过程和通路的基因，比较完全。或者是每一种里的"all"。但是你如果做免疫，需要专业性更强的database，那么就选C7。总之，根据你的实验需求来选。这里我选择H。
接下来，“number of permutation",置换次数。1000是默认值。内存够大CPU够强，就选1000，这里选100。越大越好，越大结果越精确。电脑不行的请选择100，当然了值越大，出来的结果会越精确。参考：学习基因通路富集分析软件GSEA
更多的置换次数需要更长的计算时间。 为了计算每个gene set的FDR Q值，通过置换每个基因组中的基因并重新计算随机组的P值来随机化数据集,此参数指定完成此随机化的次数。执行的随机化越多，FDR Q值估计就越精确（达到极限，因为最终FDR Q值将稳定在实际值）。
然后，phenotype label:
我选择的是tumor vs normal。
下面的collapse/Remap to gene symbols的选择，如果你的列表里基因名称和我的一样，都是ENSG开头的，那么你需要选择collapse；如果你的基因名称是gene symbol，那就选No_collapse。
Permutaion type的选择：样品数量少，选gene_set；样品数量多，选phenotype。
这一部分的最后一个参数：chip platform。点开也有很多的选择，有的文章里说这一步不用选，但是对我来说，如果我空着它，运行时会一直报错！这里我选择的是human_ensembl_gene_ID这一项：
下面是第二部分的参数，我是都用的默认值，只是修改一下输出文件的储存路径：
最后点击最下方的"Run"，如果你的文件格式没有错，选择的参数也对应你的文件，比如你的基因名称和collapse/no_collapse是否对应。那么在左下角就会出现running的字样：
等一段时间（取决于你的电脑性能，以及你选择的参数），运行成功会显示：
NOTE:如果点击运行，左下角出现红色字Error，那就是你的文件格式有问题，要么就是你collapse选择错了，要么就是你空着chip platform了。
（4）生成的文件怎么看
运行后，生成了一堆的文件：
这些文件那些比较重要呢？首先你可以找一下这个文件：
点开它，看一下summary：
然后我们点开看看详细情况，会告诉你这24个gene set都是什么：
点击details,会出现对于这个gene set的分析结果:
重点来了：这个图怎么看？
先看官网放的图，有点糊：
根据官网的图，我得到的这个图分为三部分：
（1）第一部分：绿色的折线图。就是我们上面说到的加分减分得到的曲线。我的图这里的峰是向下的，它的富集分数ES是-0.733。
（2）第二部分：中间黑色的像条形码一样的图。是这个gene set里所有的基因，这个gene set里有200个基因，所以这些黑线理论上有200条。可以看到大部分的黑线都集中在右边，也就是说这个gene set里的大部分基因在tumor vs normal里是下调的。
（3）第三部分：灰色部分。是所有基因的rank 值分布图。
你还会在下面看到这个gene set里的所有基因的情况：
关于上面的有些参数，我就不详细说明了，请参考文章，这两篇文章里有相当详细的说明，我就不做搬运工了，copy也没啥意思：
1.学习基因通路富集分析软件GSEA
2.GSEA教程
总结
学习生信一定要自己亲自实践一遍！光看是没有用的！因为网上的教程很多，说的步骤也都是五花八门，不同文章里的不同点就是你在操作的时候容易出错的地方。我也跟着不同的教程走流程，发现不是每一个教程都完全正确。在操作过程中会遇到各种各样的问题，报错报到头大。有时一个报错要各种的搜索文章，各种百度google，折腾上好几天也是常有的事。但是每次靠着自己解决问题的那种成就感是无法言喻的。所以要相信自己！只要没气到把电脑砸了，就一定会找到解决问题的方法！
加油！
