主讲老师第三十三课：GSVA分析
对于富集分析的标准流程，包括 GO 分析，KEGG 分析，GSEA 分析，然而，标准分析通常是不够的，定位到成百上千个有统计学显著变化的差异表达基因后，同样是有成百上千个生物学功能注释（GO 功能和 KEGG 通路），当普通的分析已经不能满足大家多元化的分析时，所以就有了大家耳熟能详的GSEA 分析，以及绝大部分人比较陌生的 GSVA 分析（GSVA.zip），GSVA 其实就是计算每个 sample 的GSEA 然后得出类似 pathway enrich score，把这个可以当作芯片的表达数据一样，再用 limma 包分析
差异基因。GSVA 分析的文章发表
于 2013 年，题目是 GSVA: gene set variation analysis for microarray and RNA-Seq data然而，同样是 broad 研究生出品，GSVA 文章的引用率远远不如其在 2005 年 PNAS 发表的 gsea 的引用。接下来，讲一下 GSVA 分析，一般人做数据挖掘，到差异基因的生物学功能注释（GO 功能和KEGG 通路）就结束了，进而也就是去使用一些网页工具，比如 string，出一些花花绿绿的图表，比如 PPI 网络图，实际上，使用了 GSVA，其实就是 pathway 级别的差异分析，可以把成百上千个生物学功能注释（GO 功能和 KEGG 通路）转换为新的表达矩阵，即具体的每个通路在各个样本的基因集变异分析（Gene Set Variation Analysis，GSVA） 值。下面，我们来看下如何
进行 GSVA 分析。
1.读取标准化后的表达矩阵
首先，读取准备好的 TPM 表达数据，并使用 as.matrix()函数将其转换为矩阵形式，对了，这部分代码千万别在课上运行，不然会进退两难。通过 dim()函数，可以看到，其中包含 19584 个基因表达信息，577 例样品
2.R 包与文件准备
2.1 加载 R 包
加载 GSVA 包进行后续的分析。
2.2 读取基因列表文件
接着，与 GSEA 分析类似，我们需要设置 gmt 文件路径，与 GSVA 分析不同，这次我们使用
msigdb.v7.0 数据集进行分析
然后，使用 getGmt()函数获取 gmt 文件内容。
3.GSVA 分析
准备工作完成后，接下来就可以开始正式的 GSVA 分析。
在运行过程中，会计算每个样品在每条通路中的富集分数，最终输出一个行为通路，列为样品的矩阵，因为数据量较大，这一句代码大概要运行半个小时。
最终，将 GSVA 运行结果进行保存。
4.GSVA 差异分析
4.1 分组信息
差异分析，自然首先需要设置分组信息。
在此，根据组织的来源，将其分成肿瘤组 Tumor 和正常组 Normal，并将其转换成因子形式。
结果显示：其中肿瘤组织 519 个，正常组织 58 个。
4.2 差异分析
随后，使用 limma 包进行差异分析，此处，差异分析方法与 GEO 芯片的差异分析方法类似，基于分组信息 group_list，设置比较方法为 Tumor vs. Normal。
然后，进一步拟合 limma 线性模型定义 P 值=0.05 为分界，对结果的显著性水平进行分析结果显示：其中，下调的 12598 条通路，上调的 7368 条通路。
使用 topTable()函数提取差异分析结果，并使用 na.omit()函数去除包含 NA 的行。
最后，将差异分析结果和分组信息进行保存。
5.可视化分析
首先，将分析得到的 GSVA 矩阵和差异分析结果读取进来
根据 P 值小于 0.05，且 logFC 大于 0.5，进行筛选，最终得到 10 条显著差异的富集通路。
随后，对富集得到的显著差异表达通路进行提取子集，用于绘制热图的输入文件根据分组情况，构建注释文件，使用 pheatmap()函数进行热图的绘制可以看到，显著富集的几条通路在肿瘤中基本呈现下调的结果，到此，富集分析就先介绍到这里了。