/
gsva.Rmd
182 lines (131 loc) · 8.52 KB
/
gsva.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
---
title: "GSVA: gene set variation analysis"
date: "`r Sys.Date()`"
output:
workflowr::wflow_html:
toc: true
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
Following the [vignette](https://bioconductor.org/packages/release/bioc/vignettes/GSVA/inst/doc/GSVA.html).
>Gene set variation analysis (GSVA) is a particular type of gene set enrichment method that works on single samples and enables pathway-centric analyses of molecular data by performing a conceptually simple but powerful change in the functional unit of analysis, from genes to gene sets. The GSVA package provides the implementation of four single-sample gene set enrichment methods, concretely zscore, plage, ssGSEA and its own called GSVA. While this methodology was initially developed for gene expression data, it can be applied to other types of molecular profiling data. In this vignette we illustrate how to use the GSVA package with bulk microarray and RNA-seq expression data.
>
>Gene set variation analysis (GSVA) provides an estimate of pathway activity by transforming an input gene-by-sample expression data matrix into a corresponding gene-set-by-sample expression data matrix.
>
>This resulting expression data matrix can be then used with classical analytical methods such as differential expression, classification, survival analysis, clustering or correlation analysis in a pathway-centric manner. One can also perform sample-wise comparisons between pathways and other molecular data types such as microRNA expression or binding data, copy-number variation (CNV) data or single nucleotide polymorphisms (SNPs).
## Package
Install GSVA. (Dependencies are listed in the Imports section in the [DESCRIPTION](https://github.com/rcastelo/GSVA/blob/devel/DESCRIPTION) file.)
```{r install_fgsea, message=FALSE, warning=FALSE}
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if (!require("GSVA", quietly = TRUE))
BiocManager::install("GSVA")
```
Load package.
```{r load_gsva, message=FALSE, warning=FALSE}
library(GSVA)
packageVersion("GSVA")
```
## Quick start
Generate example expression matrix.
```{r eg_mat}
p <- 10000
n <- 30
set.seed(1984)
X <- matrix(
rnorm(p*n),
nrow=p,
dimnames=list(paste0("g", 1:p), paste0("s", 1:n))
)
X[1:5, 1:5]
```
Generate 100 gene sets that are contain from 10 to up to 100 genes sampled from `1:p`.
```{r eg_gs}
set.seed(1984)
gs <- as.list(sample(10:100, size=100, replace=TRUE))
gs <- lapply(gs, function(n, p){
paste0("g", sample(1:p, size=n, replace=FALSE))
}, p)
names(gs) <- paste0("gs", 1:length(gs))
sapply(gs, length)
```
Calculate GSVA enrichment scores using the `gsva()` function, which does all the work and requires the following two input arguments:
1. A **normalised** gene expression dataset, which can be provided in one of the following containers:
* A matrix of expression values with genes corresponding to rows and samples corresponding to columns.
* An ExpressionSet object; see package Biobase.
* A SummarizedExperiment object, see package SummarizedExperiment.
2. A collection of gene sets; which can be provided in one of the following containers:
* A list object where each element corresponds to a gene set defined by a vector of gene identifiers, and the element names correspond to the names of the gene sets.
* A GeneSetCollection object; see package GSEABase.
>The first argument to the `gsva()` function is the gene expression data matrix and the second the collection of gene sets. The `gsva()` function can take the input expression data and gene sets using different specialized containers that facilitate the access and manipulation of molecular and phenotype data, as well as their associated metadata. Another advanced features include the use of on-disk and parallel backends to enable, respectively, using GSVA on large molecular data sets and speed up computing time.
The `gsva()` function will apply the following filters before the actual calculations take place:
1. Discard genes in the input expression data matrix with constant expression.
2. Discard genes in the input gene sets that do not map to a gene in the input gene expression data matrix.
3. Discard gene sets that, after applying the previous filters, do not meet a minimum and maximum size, which by default is 1 for the minimum size and `Inf` for the maximum size.
When `method="gsva"` is used (the default), the following parameters can be tuned:
* `kcdf`: The first step of the GSVA algorithm brings gene expression profiles to a common scale by calculating an expression statistic through a non-parametric estimation of the CDF across samples. Such a non-parametric estimation employs a kernel function and the `kcdf` parameter allows the user to specify three possible values for that function:
1. "Gaussian", the default value, which is suitable for continuous expression data, such as microarray fluorescent units in logarithmic scale and RNA-seq log-CPMs, log-RPKMs or log-TPMs units of expression;
2. "Poisson", which is suitable for integer counts, such as those derived from RNA-seq alignments;
3. "none", which will enforce a direct estimation of the CDF without a kernel function.
* `mx.diff`: The last step of the GSVA algorithm calculates the gene set enrichment score from two Kolmogorov-Smirnov random walk statistics. This parameter is a logical flag that allows the user to specify two possible ways to do such calculation:
1. `TRUE`, the default value, where the enrichment score is calculated as the magnitude difference between the largest positive and negative random walk deviations;
2. `FALSE`, where the enrichment score is calculated as the maximum distance of the random walk from zero.
* `abs.ranking`: Logical flag used only when `mx.diff=TRUE`. By default, `abs.ranking=FALSE` and it implies that a modified Kuiper statistic is used to calculate enrichment scores, taking the magnitude difference between the largest positive and negative random walk deviations. When `abs.ranking=TRUE` the original Kuiper statistic is used, by which the largest positive and negative random walk deviations are added together. In this case, gene sets with genes enriched on either extreme (high or low) will be regarded as highly activated.
* `tau`: Exponent defining the weight of the tail in the random walk. By default `tau=1`. When `method="ssgsea"`, this parameter is also used and its default value becomes then `tau=0.25` to match the methodology described in (Barbie et al. 2009).
In general, the default values for the previous parameters are suitable for most analysis settings, which usually consist of some kind of normalized continuous expression values.
```{r eg_gsva}
es_gsva <- gsva(gsvaParam(X, gs), verbose=FALSE)
dim(es_gsva)
```
Median enrichment scores.
```{r eg_median_es}
apply(es_gsva, 2, median)
```
>ssgsea (Barbie et al. 2009). Single sample GSEA (ssGSEA) is a non-parametric method that calculates a gene set enrichment score per sample as the normalized difference in empirical cumulative distribution functions (CDFs) of gene expression ranks inside and outside the gene set. By default, the implementation in the GSVA package follows the last step described in (Barbie et al. 2009, online methods, pg. 2) by which pathway scores are normalized, dividing them by the range of calculated values. This normalization step may be switched off using the argument ssgsea.norm in the call to the gsva() function; see below.
```{r eg_ssgsea}
es_ssgsea <- gsva(ssgseaParam(X, gs), verbose=FALSE)
apply(es_ssgsea, 2, median)
```
## Investigating the scores
Create another test matrix.
```{r create_test_matrix}
p <- 10000
n <- 2
set.seed(1984)
X <- matrix(
rnorm(n = p*n, mean = 10, sd = 10),
nrow=p,
dimnames=list(paste0("g", 1:p), paste0("s", 1:n))
)
X[1:50, 's1'] <- rnorm(n = 50, mean = 50, sd = 55)
X[51:100, 's1'] <- rnorm(n = 50, mean = 2, sd = 2)
X[1:50, 's2'] <- rnorm(n = 50, mean = 100, sd = 5)
X[51:100, 's2'] <- rnorm(n = 50, mean = 2, sd = 0.5)
X[1:5, ]
```
Create testing gene lists with higher and lower expression patterns and check their enrichment scores.
```{r ssgsea_test_matrix}
gene_set <- list(
gs1 = paste0("g", 1:50),
gs2 = paste0("g", 51:100),
gs3 = paste0("g", 101:150)
)
test_ssgsea <- gsva(ssgseaParam(X, gene_set), verbose=FALSE)
test_ssgsea
```
GSVA method.
```{r gsva_test_matrix}
test_gsva <- gsva(gsvaParam(X, gene_set), verbose=FALSE)
test_gsva
```
Z-score method.
```{r zscore_test_matrix}
test_zscore <- gsva(zscoreParam(X, gene_set), verbose=FALSE)
test_zscore
```
PLAGE method.
```{r plage_test_matrix}
test_plage <- gsva(plageParam(X, gene_set), verbose=FALSE)
test_plage
```