PGSEA和GSVA你会怎么选择呢？
GSEA 相信看过我生信菜鸟团博客的朋友都已经耳熟能详了的，其需要样本的描述以及分组信息，来计算每个基因的差异度量对它们进行排序，然后走GSEA。 虽然有ssGSEA这样的单样本的分析，但仍然不够，也有GSVA这样的算法来弥补，这里要介绍的是另外一个包，PGSEA。
it tests for each sample whether the average expression of genes in a gene sets deviates from the overall average expression (expression of all genes in all samples).
使用GSVA方法计算某基因集在各个样本的表现
安装PGSEA这个R包
安装并且查看 PDF教程：
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("PGSEA")
library(PGSEA)
browseVignettes("PGSEA") 最新版教程： https://bioconductor.org/packages/release/bioc/html/PGSEA.html
重点就是PGSEA函数及smcPlot函数，前者根据表达矩阵及基因集来进行GSEA分析，后者用来可视化分析后的结果。
因为PGSEA分析后的结果是每个基因集在每个样本的一个score，所以也是一个表达矩阵，也可以进行limma的差异分析流程。
执行PGSEA分析及可视化
分析结果也可以走limma流程
library(PGSEA)
library(GEOquery)
library(GSEABase)
gse <- getGEO("GSE7023",GSEMatrix=TRUE)
#load("gse.rda")
subtype <- gsub("\\.", "_",gsub("subtype: ", "", phenoData(gse[[1]])$"characteristics_ch1"))
pheno <- new("AnnotatedDataFrame", data = data.frame(subtype), varMetadata = data.frame(labelDescription="subtype"))
rownames(pheno@data) <- colnames(exprs(gse[[1]]))
eset <- new("ExpressionSet", exprs = exprs(gse[[1]]), phenoData = pheno)
data(VAIgsc)
details(VAIgsc[[1]])
pgNF <- PGSEA(eset, VAIgsc, ref=which(subtype=="NO"), p.value=NA)
library(limma)
design <- model.matrix(~ -1+factor(subtype))
colnames(design) <- names(table(subtype))
fit <- lmFit(pgNF, design)
contrast.matrix <- makeContrasts(P2B-NO , levels=design)
fit <- contrasts.fit(fit, contrast.matrix)
fit <- eBayes(fit)
topTable(fit, n=10)[,c("logFC","t","adj.P.Val")]可以看到都是给一个基因集（GO/KEGG/BIOCARTA/REACTOME/MSIGDB）在每个样本里面打分，把所有基因在所有样本的表达量矩阵转换为了所有基因集在所有样本的打分矩阵！！！
那么，回答我，你会怎么选择呢？
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2018-06-13，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录