一般我們做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 三種計算基因集得分方法。
作者：BeeBee生信
鏈接：http://www.reibang.com/p/1a82b63d65c3
