GSVA其实就是pathway级别的差异分析,根据GSVA的原理其实就是计算每个sample的GSEA然后得出类似pathway enrich score,把这个可以当作芯片的表达数据一样,再用limma包分析差异基因。
参考以下两个例子:
使用GSVA方法计算某基因集在各个样本的表现
充分理解GSVA和GSEA
以及我一个例子:
rm(list=ls())
suppressMessages(library(GSVA))
suppressMessages(library(GSVAdata))
suppressMessages(library(GSEABase))
suppressMessages(library(limma))
#读入gmt文件,这个可以从MSigDB上下载,这边选的上gene symbol根据自己的data来选择
gmt_file="~/c5.bp.v7.1.symbols.gmt"
geneset <- getGmt(gmt_file)
#自行读入exp文件,我这边为行为gene symbol,列为样本(4个control,5个treatment样本)
es <- gsva(as.matrix(exp), geneset,
min.sz=10, max.sz=500, verbose=TRUE)
#得到gsva计算的数值后再用limma包做差异分析得到差异的pathway
design <- model.matrix(~ factor(c(rep("cont",4),rep("treatment",5))))
colnames(es) <- c("ALL", "contvstreatment")
row.names(design)<-colnames(exp)
fit <- lmFit(es, design)
fit <- eBayes(fit)
#这边是总的
allGeneSets <- topTable(fit, coef="contvstreatment", number=Inf)
#这边是差异的
adjPvalueCutoff <- 0.001
DEgeneSets <- topTable(fit, coef="MLLvsALL", number=Inf,
p.value=adjPvalueCutoff, adjust="BH")
res <- decideTests(fit, p.value=adjPvalueCutoff)
网友评论