一般我们做GSEA都是先进行差异基因分析，然后取差异倍数排序结果进行GSEA。但如果你没有条件进行差异基因分析也可以进行GSEA，这就是ssGSEA(single sample GSEA, 单样品GSEA)对每个样品进行GSEA。R包GSVA(Gene Set Variation Analysis)可以用来进行ssGSEA。GSVA将 表达矩阵转换成通路富集分数(ES)矩阵 ，再借用limma包的 lmFit 分析得到差异通路。
先导入需要的包和检查一下表达矩阵，我例子表达矩阵是6个样品的RNA-seq数据(read count)。其中 GSEABase 包用于读取gmt格式基因集文件，本文用从MSigDB下载的KEGG基因集。
library(GSVA)
library(GSEABase)
library(limma)
> head(readCount)
A1 A2 A3 B1 B2 B3
1 807 1102 1252 237 689 485
2 50 93 28 38 2104 1083
9 292 839 319 301 192 327
10 26 8 0 0 4 0
12 3 214 41 0 1 2
13 5 0 0 9 131 75
> dim(readCount)
[1] 20292 6
因为基因集使用ENTREZID标志基因，我的表达矩阵基因名(即行名)也转换成ENTREZID，这点要保持一致。然后读入KEGG基因集，用GSVA将表达矩阵转换成通路矩阵。此时已经获得单样本的通路富集结果。这里 kcdf 参数设为"Poisson"是因为使用read count数据，如果是使用 log 后的CPM, RPKM, TPM等数据就用默认值"Gaussian"。参数 parrallel.sz 设置并行线程数，因为每个基因集的计算是独立的。
> keggSet <- getGmt("/datapool/pengguoyu/Database/MSigDB/c2.cp.kegg.v7.0.entrez.gmt")
> keggEs <- gsva(expr=as.matrix(readCount), gset.idx.list=keggSet, kcdf="Poisson", parallel.sz=4)
> head(keggEs, n=3)
A1 A2 A3 B1
KEGG_GLYCOLYSIS_GLUCONEOGENESIS -0.1861730 0.2241277 0.0004014829 -0.3203273
KEGG_CITRATE_CYCLE_TCA_CYCLE -0.5081698 0.5308383 0.1455315128 -0.4934121
KEGG_PENTOSE_PHOSPHATE_PATHWAY -0.2435699 0.3432502 -0.2706873733 -0.1581674
B2 B3
KEGG_GLYCOLYSIS_GLUCONEOGENESIS -0.10393707 0.02900289
KEGG_CITRATE_CYCLE_TCA_CYCLE -0.25016781 0.02487492
KEGG_PENTOSE_PHOSPHATE_PATHWAY 0.07934649 -0.28966916
> dim(keggEs)
[1] 186 6
准备用limma的 lmFit 设定组间比较为B组比A组。
> grouP <- c(rep("A", 3), rep("B", 3)) %>% as.factor()
> desigN <- model.matrix(~ grouP + 0)
> rownames(desigN) <- c("A1", "A2", "A3", "B1", "B2", "B3")
> desigN
grouPA grouPB
A1 1 0
A2 1 0
A3 1 0
B1 0 1
B2 0 1
B3 0 1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$grouP
[1] "contr.treatment"
# 用desigN的列名，设定B组比A组
> comparE <- makeContrasts(grouPB - grouPA, levels=desigN)
使用limma取得差异通路结果，这里KEGG一共才186通路，所以 number 参数设200就可以取得所有结果。
> fiT <- lmFit(keggEs, desigN)
> fiT2 <- contrasts.fit(fiT, comparE)
> fiT3 <- eBayes(fiT2)
> keggDiff <- topTable(fiT3, coef=1, number=200)
> head(keggDiff, n=3)
logFC AveExpr
KEGG_GLYCOSAMINOGLYCAN_DEGRADATION -0.7343229 -0.006312659
KEGG_FOCAL_ADHESION -0.6204141 -0.056091519
KEGG_GLYCOSAMINOGLYCAN_BIOSYNTHESIS_CHONDROITIN_SULFATE -0.8470823 -0.039176682
t P.Value
KEGG_GLYCOSAMINOGLYCAN_DEGRADATION -6.643722 8.986989e-05
KEGG_FOCAL_ADHESION -5.918378 2.147451e-04
KEGG_GLYCOSAMINOGLYCAN_BIOSYNTHESIS_CHONDROITIN_SULFATE -5.795434 2.505992e-04
adj.P.Val B
KEGG_GLYCOSAMINOGLYCAN_DEGRADATION 0.01415655 1.8704201
KEGG_FOCAL_ADHESION 0.01415655 1.0675346
KEGG_GLYCOSAMINOGLYCAN_BIOSYNTHESIS_CHONDROITIN_SULFATE 0.01415655 0.9232913
这个GSVA还支持 ssGSEA, zscore, plage 三种计算基因集得分方法。
