使用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")

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-05-29

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏iOSDevLog

Turi Create 机器学习模型实战:你也能轻松做出Prisma 风格的图片!

如果你一直有关注Apple去年所发布的消息,就会知道他们在机器学习上投入了大量心力。自他们去年在WWDC 2017上推出Core ML以来,已经有大量结合机器学...

39620
来自专栏数据小魔方

带实际执行进度的甘特图

今天要跟大家分享的图标是带实际执行进度的甘特图! ▽▼▽ 由于本图所用到的技巧和思路特别复杂,过程相对繁琐,所以本案例的介绍会省略掉很多细节性的步骤,否则图文会...

38950
来自专栏一心无二用,本人只专注于基础图像算法的实现与优化。

SSE图像算法优化系列二:高斯模糊算法的全面优化过程分享(二)。

      相关链接: 高斯模糊算法的全面优化过程分享(一)      在高斯模糊算法的全面优化过程分享(一)一文中我们已经给出了一种相当高性能的高斯模糊过程,...

45360
来自专栏PPV课数据科学社区

数据挖掘系列(4)使用weka做关联规则挖掘

前面几篇介绍了关联规则的一些基本概念和两个基本算法,但实际在商业应用中,写算法反而比较少,理解数据,把握数据,利用工具才是重要的,前面的基础篇是对算法的理解,这...

52060
来自专栏章鱼的慢慢技术路

用OpenGL实现粒子的随机运动

17150
来自专栏生信技能树

【直播】我的基因组51:画全基因范围内的染色体reads覆盖度图

前面我们已经详细讲解过如何根据窗口来统计每条染色体的每个片段的GC含量,还有平均测序深度,请大家自行前往前面查看脚本及实现方式!【直播】我的基因组47:测序深度...

53060
来自专栏ATYUN订阅号

【学术】不懂神经网络?不怕,一文教你用JavaScript构建神经网络

AiTechYun 编辑:xiaoshan.xiang 本文的内容并不是关于神经网络的深度教程,在这里既不会深入研究输入层、激活函数的内部原理,也不会教你如何使...

35240
来自专栏生信宝典

R包reshape2,轻松实现长、宽数据表格转换

本文翻译自外文博客,原文链接:https://seananderson.ca/2013/10/19/reshape/

25520
来自专栏Y大宽

ToppGene Suite中文使用指南

2007.12:Improved human disease candidate gene prioritization using mouse phenoty...

20230
来自专栏UAI人工智能

[译] TensorFlow 白皮书

53760

扫码关注云+社区

领取腾讯云代金券