使用GSVA方法计算某基因集在各个样本的表现

文章发表于2013年，GSVA: gene set variation analysis for microarray and RNA-Seq data 同样是broad 研究生出品，其在2005年PNAS发表的gsea已经高达1.4万的引用了，不过这个GSVA才不到300。

算法细节

算法本身就不是很好理解，并不强求一定要理解透彻，可以参考2005年的GSEA算法：

GSEA 算法

GSEA分析一文就够（单机版+R语言版）

GSEA的统计学原理试讲

GSVA starts by evaluating whether a gene i is highly or lowly expressed in sample j in the context of the sample population distribution.

可以是芯片杂交的信号代表的表达量，也可以是转录组测序定量。

For each gene expression profile xi={xi1,…,xin}, a non-parametric kernel estimation of its cumulative density function is performed.

We offer two approaches for turning the KS like random walk statistic into an enrichment statistic (ES) (also called GSVA score), the classical maximum deviation method and a normalized ES.

而且作者也在测试数据和真实数据把自己的GSVA算法跟GSEA，PLAGE, single sample GSEA (ssGSEA)或者其它算法进行了比较， 还在TCGA的ovarian serous cystadenocarcinoma (OV)癌症表达矩阵(n=588) ，用MSigDB数据库的 canonical gene sets (C2) 基因集做了比较和测试。

还比较了转录组测序数据和芯片数据，这些数据都提供了下载链接，最后作者把算法打包成了 Bioconductor package for R under the name GSVA at http://www.bioconductor.org.

安装GSVA这个R包

安装并且查看21页的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("GSVA")
library(GSVA)
browseVignettes("GSVA")
browseVignettes("estimate")


最新版教程：https://www.bioconductor.org/packages/devel/bioc/vignettes/GSVA/inst/doc/GSVA.pdf

其实核心函数就是gsva()，需要两个输入：the gene expression data and a collection of gene sets.

其实这个函数也可以选择其它3个模型：

• method="plage" (Tomfohr et al., 2005). Pathway level analysis of gene expression (PLAGE)

• method="zscore" (Lee et al., 2008). The combined z-score method a

• method="ssgsea" (Barbie et al., 2009). Single sample GSEA (ssGSEA) calculates a gene set
  enrichment score per sample

另外一个比较重要的参数是： default argument mx.diff=TRUE to obtain approximately
normally distributed ES，如果设置为false，那么通常是 a bimodal distribution of GSVA enrichment scores for each gene

非常多的文章都在引用该算法，比如：https://www.nature.com/articles/srep16238#f1

先在模拟数据应用GSVA

代码很简单，构造一个 30个样本，2万个基因的表达矩阵， 加上 100 个假定的基因集。

library(GSVA)

p <- 20000    ## number of genes
n <- 30       ## number of samples
nGS <- 100    ## number of gene sets
min.sz <- 10  ## minimum gene set size
max.sz <- 100 ## maximum gene set size
X <- matrix(rnorm(p*n), nrow=p, dimnames=list(1:p, 1:n))
dim(X)
gs <- as.list(sample(min.sz:max.sz, size=nGS, replace=TRUE)) ## sample gene set sizes
gs <- lapply(gs, function(n, p) sample(1:p, size=n, replace=FALSE), p) ## sample gene sets
es.max <- gsva(X, gs, mx.diff=FALSE, verbose=FALSE, parallel.sz=1)
es.dif <- gsva(X, gs, mx.diff=TRUE, verbose=FALSE, parallel.sz=1)
pheatmap::pheatmap(es.max)
pheatmap::pheatmap(es.dif)


这样就可以检验我们假定的100个基因集在我们的30个样本的GSVA score值分布情况。

值得注意的是，这里的gsva函数接受的是一个纯粹的表达矩阵matrix和一个纯粹基因集合list，实际上通常是一个 ExpressionSet 和 GeneSetCollection 对象，所以大家务必学会R里面的对象哈，不然永远只能看懂简单代码。

两个模型的区别

par(mfrow=c(1,2), mar=c(4, 4, 4, 1))
plot(density(as.vector(es.max)), main="Maximum deviation from zero",
    xlab="GSVA score", lwd=2, las=1, xaxt="n", xlim=c(-0.75, 0.75), cex.axis=0.8)
axis(1, at=seq(-0.75, 0.75, by=0.25), labels=seq(-0.75, 0.75, by=0.25), cex.axis=0.8)
plot(density(as.vector(es.dif)), main="Difference between largest\npositive and negative deviations",
    xlab="GSVA score", lwd=2, las=1, xaxt="n", xlim=c(-0.75, 0.75), cex.axis=0.8)
axis(1, at=seq(-0.75, 0.75, by=0.25), labels=seq(-0.75, 0.75, by=0.25), cex.axis=0.8)
# 前者是高斯分布，后者是二项式分布


在真实数据(白血病芯片数据)上面测试GSVA算法

library(GSEABase)
library(GSVAdata)
data(c2BroadSets)
c2BroadSets
head(names(c2BroadSets))
c2BroadSets[[names(c2BroadSets)[1]]]

library(Biobase)
library(genefilter)
library(limma)
library(RColorBrewer)
library(GSVA)

cacheDir <- system.file("extdata", package="GSVA")
cachePrefix <- "cache4vignette_"
file.remove(paste(cacheDir,
                 list.files(cacheDir, pattern=cachePrefix), sep="/"))

data(leukemia)
leukemia_eset
head(pData(leukemia_eset))
table(leukemia_eset$subtype)


是17个MLL病人和20个ALL病人的古老的affymetrix芯片表达数据，共12626个探针，这样的表达矩阵首先需要过滤哈，这里直接用 genefilter 包提供的 nsFilter函数来进行过滤。

filtered_eset <- nsFilter(leukemia_eset, require.entrez=TRUE,
                         remove.dupEntrez=TRUE,
                         var.func=IQR, var.filter=TRUE,
                         var.cutoff=0.5, filterByQuantile=TRUE,
                         feature.exclude="^AFFX")
filtered_eset
leukemia_filtered_eset <- filtered_eset$eset
leukemia_filtered_eset
exprs(leukemia_filtered_eset)[1:4,1:4]


剩下的是 4292个探针的表达矩阵。

根据表型数据使用limma包来找到有显著差异的基因集

因为每个基因集都在每个样本里面得到了一个值，所以这时候相当于有了一个新的表达矩阵，而且这些样本的表型数据仍然是存在的，所以可以借鉴差异分析的算法了。

cache(leukemia_es <- gsva(leukemia_filtered_eset, c2BroadSets,
                         min.sz=10, max.sz=500, verbose=TRUE),
     dir=cacheDir, prefix=cachePrefix)

adjPvalueCutoff <- 0.001
logFCcutoff <- log2(2)
design <- model.matrix(~ factor(leukemia_es$subtype))
colnames(design) <- c("ALL", "MLLvsALL")
fit <- lmFit(leukemia_es, design)
fit <- eBayes(fit)
allGeneSets <- topTable(fit, coef="MLLvsALL", number=Inf)
DEgeneSets <- topTable(fit, coef="MLLvsALL", number=Inf,
                      p.value=adjPvalueCutoff, adjust="BH")
res <- decideTests(fit, p.value=adjPvalueCutoff)
summary(res)


Thus, there are 38 MSigDB C2 curated pathways that are differentially activated between MLL and ALL
at 0.1% FDR.

同样limma当然是可以找显著差异表达的基因的

logFCcutoff <- log2(2)
design <- model.matrix(~ factor(leukemia_eset$subtype))
colnames(design) <- c("ALL", "MLLvsALL")
fit <- lmFit(leukemia_filtered_eset, design)
fit <- eBayes(fit)
allGenes <- topTable(fit, coef="MLLvsALL", number=Inf)
DEgenes <- topTable(fit, coef="MLLvsALL", number=Inf,
                   p.value=adjPvalueCutoff, adjust="BH", lfc=logFCcutoff)
res <- decideTests(fit, p.value=adjPvalueCutoff, lfc=logFCcutoff)
summary(res)


Here, 122 genes show up as being differentially expressed with a minimum fold-change of 2 at 0.1% FDR.

可以看到，两个代码唯一的变化就是 leukemia_filtered_eset 和 leukemia_es而已。这样的差异分析结果同样也是可以画火山图，热图，代码就不赘述了，非常简单。

先看两个火山图的区别：

然后看两个热图的区别；

还可以对包内置的TCGA数据进行测试

data(gbm_VerhaakEtAl)
gbm_eset
exprs(gbm_eset)[1:4,1:4]
### 如下
ExpressionSet (storageMode: lockedEnvironment)
assayData: 11861 features, 173 samples
 element names: exprs
protocolData: none
phenoData
 rowNames: TCGA.02.0003.01A.01 TCGA.02.0010.01A.01
   ... TCGA.12.0620.01A.01 (173 total)
 varLabels: subtype
 varMetadata: labelDescription channel
featureData: none
experimentData: use 'experimentData(object)'
Annotation:  
>


同样的方式做gsva分析

data(gbm_VerhaakEtAl)
gbm_eset
head(featureNames(gbm_eset))
table(gbm_eset$subtype)
data(brainTxDbSets)
sapply(brainTxDbSets, length)
lapply(brainTxDbSets, head)
## 可以看到，都是symbol格式的基因ID
gbm_es <- gsva(gbm_eset, brainTxDbSets, mx.diff=FALSE, verbose=FALSE, parallel.sz=1)


不同算法在转录组测序数据的表现

前面我们说到过gsva函数还提供了另外3个算法，这里就不细细讲解了。

也可以只关注3个主流通路数据库：

canonicalC2BroadSets <- c2BroadSets[c(grep("^KEGG", names(c2BroadSets)),
                                    grep("^REACTOME", names(c2BroadSets)),
                                     grep("^BIOCARTA", names(c2BroadSets)))]
canonicalC2BroadSets


还可以增加自己感兴趣的基因集

data(genderGenesEntrez)

MSY <- GeneSet(msYgenesEntrez, geneIdType=EntrezIdentifier(),
              collectionType=BroadCollection(category="c2"), setName="MSY")
MSY
XiE <- GeneSet(XiEgenesEntrez, geneIdType=EntrezIdentifier(),
              collectionType=BroadCollection(category="c2"), setName="XiE")
XiE

canonicalC2BroadSets <- GeneSetCollection(c(canonicalC2BroadSets, MSY, XiE))
canonicalC2BroadSets


也可以去broad官网里面下载最新版gmt文件，然后读入到R语言

myC7 <- getGmt("c7.all.v5.1.entrez.gmt", geneIdType=EntrezIdentifier(), collectionType=BroadCollection(category="c7"), sep="\t")


这个时候需要详细了解 GSEABase 包的设计。

library(GSEABase)
library(GSVAdata)
data(c2BroadSets)
c2BroadSets
head(names(c2BroadSets))
c2BroadSets[[names(c2BroadSets)[1]]]
geneIds(c2BroadSets[[names(c2BroadSets)[1]]])
### 假设之前是其它ID
library(hgu95av2.db)
getEG(geneIds(c2BroadSets[[names(c2BroadSets)[1]]]),"hgu95av2")