单样本基因集富集分析——ssGSEA

1. 概述

单样本基因集富集分析（single sample gene set enrichment analysis, ssGSEA），是GSEA方法的扩展，计算每个样本和基因集配对的富集分数。每个ssGSEA富集评分代表了样本中特定基因集的成员被协调上调或下调的程度。相关文章于2009年发表于Nature：

Title: Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1。

Year: 2009

Journal: Nature

DOI: 10.1038/nature08460

GSEA生成一个基因集的富集得分，与一个数据集中的样本集合的表型差异有关。ssGSEA不同于以组别为单位（比如cancer vs normal）的GSEA分析，而是每个样本都可以得到相应基因集的评分。通过这种方式，ssGSEA将单个样本的基因表达谱转换为基因集富集谱。这种转化使研究人员能够根据生物过程和途径的活性水平而不是通过单个基因的表达水平来描述细胞状态。因此，ssGESA如果使用的是免疫细胞marker相关的基因集，就可以计算免疫细胞浸润评分。

具体描述可以参考：

​ （1）ssGSEAProjection

​ （2）基因或蛋白富集分析，gsea与ssgsea

2. 原理算法

算法概括来说就是首先对给定样本的基因表达值进行秩次标准化，然后利用经验累积分布函数计算富集分数（ES）。

​ a. 首先假设有一个样本的表达数据，那么他应该是这样的，第一列为基因，第二列为表达值，这样的两列的数据矩阵。首先对样本的所有基因的表达水平进行排序获得其在所有基因中的秩次rank，这些基因的集合为BG。

​ b. 假设要对其进行KEGG的分析，首先我们需要在GSEA官网找到KEGG对应的gmt文件。gmt文件主要格式是：每行表示一个通路，第一列为通路ID，第二列为通路对应的描述，第三列开始到最后一列为该通路中的基因。

​ c. 那么对于任意的一个通路A，我们可以拿到这个通路的基因列表GL，从GL中寻找BG里存在的基因并计数为NC，并将这些基因的表达水平加和为SG。

开始计算ES：对于任意一个表达谱中的基因 G，如果G是集合GL中的基因则他的ES等于该基因的表达水平除以SG，否则记该基因的ES等于1除以（基因集合BG总个数减去NC）

依次计算每个BG中的基因的ES值，找到其中绝对值最大的ES作为通路A的A.ES

​ d. 到此通路A的ES计算完毕，需要一个统计学方法来评估该ES是否是显著的，即非随机的。

按照上述计算ES的方法，先随机打乱表达谱中基因的表达顺序，然后再依次计算ES值，如此重复一千次，得到一千个ES值，我们根据这一千个ES值的分布，来计算A.ES在这个分布中所处的位置及出现在该位置时的概率即得到了p值 。依次分别计算每个通路的ES及p值，然后使用多重检验矫正得到每个通路的FDR。

3. 实现代码

  • 基因集需要是list为对象。
  • 基因的表达量需要是矩阵，行为基因，列为样本
  • 默认情况下，kcdf=”Gaussian”，适用于输入表达式值连续的情况，如对数尺度的微阵列荧光单元、RNA-seq log-CPMs、log-RPKMs或log-TPMs。当输入表达式值是整数计数时，比如那些从RNA-seq实验中得到的值，那么这个参数应该设置为kcdf=”Poisson”
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
library(GSEABase)
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GSVA")
library(GSVA)

#1. 获取geneSets
# 下载28种免疫细胞的参考基因集 <http://cis.hku.hk/TISIDB/data/download/CellReports.txt>
geneSet <- read.csv("CellReports.txt",header = F,sep = "\t",) # 用EXCEL打开删除NA列
class(geneSet)
geneSet <- geneSet %>%
  column_to_rownames("V1")%>%t()
a <- geneSet
a <- a[1:nrow(a),]
set <- colnames(a)
l <- list()
#i <- "Activated CD8 T cell"
for (i in set) {
  x <-  as.character(a[,i])
  x <- x[nchar(x)!=0]
  x <-  as.character(x)
  l[[i]] <-x
}
save(l,file = "./gene_set.Rdata")

#2. 读取基因的表达矩阵
rm(list=ls())
#library(limma)
## 载入数据
load(file = "./exprSet_tpm.Rdata")#这是之前处理的一个LIHC的tpm数据
exprSet_tpm <- exprSet_tpm[,-1]
exprSet_tpm <- exprSet_tpm[,c(-425,-426,-427)]
logTPM <- log(exprSet_tpm)
load(file = "./gene_set.Rdata")

#3. 开始进行ssGSEA
library(GSVA)
ssgsea<- gsva(dat, l, method='ssgsea',kcdf='Gaussian',abs.ranking=TRUE)

# Min-Max标准化是指对原始数据进行线性变换，将值映射到[0，1]之间
# 这里是将每个样本中不同的免疫细胞比例标准化到0-1之间
ssgsea.1 <- ssgsea
for (i in colnames(ssgsea)) {
  #i <- colnames(ssgsea)[1]
  ssgsea.1[,i] <- (ssgsea[,i] -min(ssgsea[,i]))/(max(ssgsea[,i] )-min(ssgsea[,i] ))
  
}
apply(ssgsea.1[,1:6], 2, range)

4. 可视化呈现

1
2
3
4
5
6
#热图
library(pheatmap)
library(RColorBrewer)
png(file="heatmap.png")
pheatmap((ssgsea.1)
dev.off()

未完待续……