一般我們做GSEA都是先進(jìn)行差異基因分析吕嘀，然后取差異倍數(shù)排序結(jié)果進(jìn)行GSEA哄酝。但如果你沒有條件進(jìn)行差異基因分析也可以進(jìn)行GSEA，這就是ssGSEA(single sample GSEA, 單樣品GSEA)對每個樣品進(jìn)行GSEA鲸匿。R包GSVA(Gene Set Variation Analysis)可以用來進(jìn)行ssGSEA弦疮。GSVA將 表達(dá)矩陣轉(zhuǎn)換成通路富集分?jǐn)?shù)(ES)矩陣 ，再借用limma包的 lmFit 分析得到差異通路寒亥。
先導(dǎo)入需要的包和檢查一下表達(dá)矩陣邮府，我例子表達(dá)矩陣是6個樣品的RNA-seq數(shù)據(jù)(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標(biāo)志基因，我的表達(dá)矩陣基因名(即行名)也轉(zhuǎn)換成ENTREZID加勤，這點要保持一致仙辟。然后讀入KEGG基因集，用GSVA將表達(dá)矩陣轉(zhuǎn)換成通路矩陣胸竞。此時已經(jīng)獲得單樣本的通路富集結(jié)果欺嗤。這里 kcdf 參數(shù)設(shè)為"Poisson"是因為使用read count數(shù)據(jù)参萄，如果是使用 log 后的CPM, RPKM, TPM等數(shù)據(jù)就用默認(rèn)值"Gaussian"卫枝。參數(shù) parrallel.sz 設(shè)置并行線程數(shù)，因為每個基因集的計算是獨立的讹挎。
> 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
準(zhǔn)備用limma的 lmFit 設(shè)定組間比較為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的列名吆玖，設(shè)定B組比A組
> comparE <- makeContrasts(grouPB - grouPA, levels=desigN)
使用limma取得差異通路結(jié)果，這里KEGG一共才186通路马篮，所以 number 參數(shù)設(shè)200就可以取得所有結(jié)果沾乘。
> 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 三種計算基因集得分方法。
