1. 前言

在生信分析过程中，我们可能会使用到测序/芯片数据，然后有一群感兴趣的基因集，现在，我们需要探讨这群基因集在不同样本中的活化情况，同时呢，我们也可以了解到，不同类型样本之间哪群样本活化/抑制了。

2. 什么是 GSVA？

官网这么说的： 基因集变异分析 (GSVA) 通过将输入的逐个基因表达数据矩阵转换为相应的逐个基因集表达数据矩阵来提供通路活动的估计。然后可以将得到的表达数据矩阵与经典分析方法一起使用，例如以通路为中心的方式进行差异表达、分类、生存分析、聚类或相关分析。还可以在通路和其他分子数据类型（例如 microRNA 表达或结合数据、拷贝数变异 (CNV) 数据或单核苷酸多态性 (SNP)）之间进行样本比较。

简单来说 就是找样本的基因和感兴趣的基因集的相关性

再说说容易弄混的 GSEA 吧 网上介绍很多，但是不够凝练，其实就一句话，GSVA 与 GSEA 不同，GSEA 先差异，再富集，GSVA 先富集，再差异（当然，学会 GSVA 之后，有时候 GSVA 也可以不做差异），这就是一个先验与后验的区别。

3. 我需要准备什么文件

当然，所需要的文件都放在公众号后台，回复相应关键字即可获取

1. 基因集，比如我这里有 4 个基因集，每个基因集中包含不同个数的基因
图片
1. 表达谱数据（需要区分标准化后的/未标准化的，芯片的/测序的），行是基因，列是样本名

 > input_matrix[1:5,1:3]
        TCGA.02.0003.01A.01 TCGA.02.0010.01A.01 TCGA.02.0011.01B.01
AACS                0.51451            -0.36739             0.27668
FSTL1              -0.28438            -1.52133            -1.02120
ELMO2               0.09697             0.32328             0.59386
CREB3L1             0.27977            -0.64587             0.37805
RPS11               0.11334             0.01263            -0.22468`</pre>

1. 样本注释文件，sample_name 顺序需要和步骤 2中一致，第一列样本名，第二列样本类型

> ann[1:5,]
          sample_name   subtype
1 TCGA.02.0003.01A.01 Proneural
2 TCGA.02.0010.01A.01 Proneural
3 TCGA.02.0011.01B.01 Proneural
4 TCGA.02.0014.01A.01 Proneural
5 TCGA.02.0024.01B.01 Proneural`</pre>

4. 开始分析

代码主要参考 GSVA 官网教程^[1] ，有任何与官网不同的地方以官网为准。

那么为啥有了官网还要自己写一篇出来呢？

官网的太全了，很多东西其实用不上
太全了就会太复杂，没有重点
有的部分可以优化，比如作图
适合自己的，才是最好的
官网的输入文件格式，感兴趣可以去观摩一下，分析是简便了，但是格式复杂化了，反而不利于实际分析

所需代码

根据输入的表达谱矩阵不同，gsva 函数中kcdf的设置也不同，简单来说：

1. RNA-seq 的原始计数用Poisson
1. 芯片微阵列或标准化后的 RNA-seq用Gaussian

当然，也可以限制基因集的大小，"min.sz=5,max.sz=500"，将其加入到gsva函数中，可以排除基因过少或过多的基因集

> library(GSVA)
library(ComplexHeatmap)
library(tidyverse)
library(circlize)

# input
input_matrix <- as.matrix(read.table("input_matrix.txt",header = T,row.names = 1,sep = "\t"))
ann = read.table("sample_ann.txt",sep = "\t",header = T)
geneset = read.table("geneSet.txt",header = F,sep = "\t",fill = T,row.names = 1)

# dataframe2list
geneset <- setNames(split(geneset, seq(nrow(geneset))), rownames(geneset))
df2list = function(x){
  x<-as.character(x)
  x <- x[x!=""]
}
geneset = lapply(geneset, df2list)

# gsva analysis
gbm_es <- gsva(gset.idx.list = geneset, expr = input_matrix, mx.diff=FALSE,
               kcdf = "Gaussian")

# result plot
column_ha = HeatmapAnnotation(foo1 = ann$subtype,
                              col = list(foo1 = c("Proneural" = "#009491",
                                                 "Neural" = "#38b282",
                                                 "Classical" = "#9bcd83",
                                                 "Mesenchymal" = "#e9e29c")))
color_set = colorRamp2(c(-1, 0, 1), c("red", "white", "blue"))
Heatmap(gbm_es,show_row_names = T, top_annotation = column_ha,
        show_column_names = F,cluster_columns = F,col = color_set)`</pre>

5. 结果

最终即可得到一张崭新的热图，横轴为每个样本，纵轴为感兴趣的基因集，上方代表不同的疾病分型（事先已知的，其实也就是sample_ann.txt中的内容，当然没有sample_ann.txt也能做GSVA，只是做出来只能看见样本间不同基因集的打分，没有表型关联），左方是算法聚类，颜色，其实也就能够看得，不同分型之间不同通路的得分情况