%\VignetteIndexEntry{Gene Set Variation Analysis} %\VignetteDepends(Biobase, methods, gplots, glmnet} %\VignetteKeywords{GSVA, GSEA, Expression, Microarray, Pathway} %\VignettePackage{GSVA} \documentclass[a4paper]{article} \usepackage{url} \usepackage{graphicx} \usepackage{longtable} \usepackage{hyperref} \usepackage{natbib} \usepackage{fullpage} \newcommand{\Rfunction}[1]{{\texttt{#1}}} \newcommand{\Robject}[1]{{\texttt{#1}}} \newcommand{\Rpackage}[1]{{\textsf{#1}}} \newcommand{\Rclass}[1]{{\textit{#1}}} \title{GSVA: The Gene Set Variation Analysis package \\ for microarray and RNA-seq data} \author{Sonja H\"anzelmann$^1$, Robert Castelo$^1$ and Justin Guinney$^2$} \begin{document} \SweaveOpts{eps=FALSE} \maketitle \begin{quote} {\scriptsize 1. Research Program on Biomedical Informatics (GRIB), Hospital del Mar Research Institute (IMIM) and Universitat Pompeu Fabra, Parc de Recerca Biom\`edica de Barcelona, Doctor Aiguader 88, 08003 Barcelona, Catalonia, Spain 2. Sage Bionetworks, 1100 Fairview Ave N., Seattle, Washington, 98109 USA } \end{quote} \begin{abstract} The \Rpackage{GSVA} package implements a non-parametric unsupervised method, called Gene Set Variation Analysis (GSVA), for assessing gene set enrichment (GSE) in gene expression microarray and RNA-seq data. In contrast to most GSE methods, GSVA performs a change in coordinate systems, transforming the data from a gene by sample matrix to a gene set by sample matrix. Thereby allowing for the evaluation of pathway enrichment for each sample. This transformation is done without the use of a phenotype, thus facilitating very powerful and open-ended analyses in a now pathway centric manner. In this vignette we illustrate how to use the \Rpackage{GSVA} package to perform some of these analyses using published microarray and RNA-seq data already pre-processed and stored in the companion experimental data package \Rpackage{GSVAdata}. \end{abstract} <>= options(width=60) pdf.options(useDingbats=FALSE) @ \section{Introduction} Gene set enrichment analysis (GSEA) \citep[see][]{mootha_pgc_1alpha_responsive_2003, subramanian_gene_2005} is a method designed to assess the concerted behavior of functionally related genes forming a set, between two well-defined groups of samples. Because it does not rely on a ``gene list'' of interest but on the entire ranking of genes, GSEA has been shown to provide greater sensitivity to find gene expression changes of small magnitude that operate coordinately in specific sets of functionally related genes. However, due to the reduced costs in genome-wide gene-expression assays, data is being produced under more complex experimental designs that involve multiple RNA sources enriched with a wide spectrum of phenotypic and/or clinical information. The Cancer Genome Atlas (TCGA) project (see \url{http://cancergenome.nih.gov}) and the data deposited on it constitute a canonical example of this situation. To facilitate the functional enrichment analysis of this kind of data, we developed Gene Set Variation Analysis (GSVA) which allows the assessment of the underlying pathway activity variation by transforming the gene by sample matrix into a gene set by sample matrix without the \textit{a priori} knowledge of the experimental design. The method is both non-parametric and unsupervised, and bypasses the conventional approach of explicitly modeling phenotypes within enrichment scoring algorithms. Focus is therefore placed on the \textit{relative} enrichment of pathways across the sample space rather than the \textit{absolute} enrichment with respect to a phenotype. The value of this approach is that it permits the use of traditional analytical methods such as classification, survival analysis, clustering, and correlation analysis in a pathway focused manner. It also facilitates sample-wise comparisons between pathways and other complex data types such as microRNA expression or binding data, copy-number variation (CNV) data, or single nucleotide polymorphisms (SNPs). However, for case-control experiments, or data with a moderate to small sample size ($<30$), other GSE methods that explicitly include the phenotype in their model are more likely to provide greater statistical power to detect functional enrichment. In the rest of this vignette we describe briefly the methodology behind GSVA, give an overview of the functions implemented in the package and show a few applications. The interested reader is referred to \citep{haenzelmann_castelo_guinney_2013} for more comprehensive explanations and more complete data analysis examples with GSVA, as well as for citing GSVA if you use it in your own work. \section{GSVA enrichment scores} A schematic overview of the GSVA method is provided in Figure \ref{methods}, which shows the two main required inputs: a matrix $X=\{x_{ij}\}_{p\times n}$ of normalized expression values (see Methods for details on the pre-processing steps) for $p$ genes by $n$ samples, where typically $p\gg n$, and a collection of gene sets $\Gamma = \{\gamma_1, \dots, \gamma_m\}$. We shall denote by $x_i$ the expression profile of the $i$-th gene, by $x_{ij}$ the specific expression value of the $i$-th gene in the $j$-th sample, and by $\gamma_k$ the subset of row indices in $X$ such that $\gamma_k \subset \{1,\ldots\,p\}$ defines a set of genes forming a pathway or some other functional unit. Let $|\gamma_k |$ be the number of genes in $\gamma_k$. \begin{figure}[ht] \centerline{\includegraphics[width=\textwidth]{methods}} \caption{{\bf GSVA methods outline.} The input for the GSVA algorithm are a gene expression matrix in the form of log2 microarray expression values or RNA-seq counts and a database of gene sets. 1. Kernel estimation of the cumulative density function (kcdf). The two plots show two simulated expression profiles mimicking 6 samples from microarray and RNA-seq. The $x$-axis corresponds to expression values where each gene is lowly expressed in the four samples with lower values and highly expressed in the other two. The scale of the kcdf is on the left $y$-axis and the scale of the Gaussian and Poisson kernels is on the right $y$-axis. 2. The expression-level statistic is rank ordered for each sample. 3. For every gene set, the Kolmogorov-Smirnov-like rank statistic is calculated. The plot illustrates a gene set consisting of 3 genes out of a total number of 10 with the sample-wise calculation of genes inside and outside of the gene set. 4. The GSVA enrichment score is either the difference between the two sums or the maximum deviation from zero. The two plots show two simulations of the resulting scores under the null hypothesis of no gene expression change (see main text). The output of the algorithm is matrix containing pathway enrichment profiles for each gene set and sample. } \label{methods} \end{figure} GSVA starts by evaluating whether a gene $i$ is highly or lowly expressed in sample $j$ in the context of the sample population distribution. Probe effects can alter hybridization intensities in microarray data such that expression values can greatly differ between two non-expressed genes\cite{zilliox_gene_2007}. Analogous gene-specific biases, such as GC content or gene length have been described in RNA-seq data\cite{hansen_removing_2012}. To bring distinct expression profiles to a common scale, an expression-level statistic is calculated as follows. For each gene expression profile $x_i=\{x_{i1},\dots,x_{in}\}$, a non-parametric kernel estimation of its cumulative density function is performed using a Gaussian kernel \cite[pg.~148]{silverman_density_1986} in the case of microarray data: \begin{equation} \label{density} \hat{F}_{h_i}(x_{ij})=\frac{1}{n}\sum_{k=1}^n\int_{-\infty}^{\frac{x_{ij}-x_{ik}}{h_i}}\frac{1}{\sqrt{2\pi}}e^{-\frac{t^2}{2}}dt\,, \end{equation} where $h_i$ is the gene-specific bandwidth parameter that controls the resolution of the kernel estimation, which is set to $h_i=s_i/4$, where $s_i$ is the sample standard deviation of the $i$-th gene (Figure \ref{methods}, step 1). In the case of RNA-seq data, a discrete Poisson kernel \cite{canale_bayesian_2011} is employed: \begin{equation} \hat{F}_r(x_{ij}) = \frac{1}{n} \sum_{k=1}^n \sum_{y=0}^{x_{ij}} \frac{e^{-(x_{ik}+r)}(x_{ik}+r)^y}{y!}\,, \end{equation} where $r=0.5$ in order to set the mode of the Poisson kernel at each $x_{ik}$, because the mode of a Poisson distribution with an integer mean $\lambda$ occurs at $\lambda$ and $\lambda-1$ and at the largest integer smaller than $\lambda$ when $\lambda$ is continuous. Let $z_{ij}$ denote the previous expression-level statistic $\hat{F}_{h_i}(x_{ij})$, or $\hat{F}_r(x_{ij})$, depending on whether $x_{ij}$ are continuous microarray, or discrete count RNA-seq values, respectively. The following step condenses expression-level statistics into gene sets by calculating sample-wise enrichment scores. To reduce the influence of potential outliers, we first convert $z_{ij}$ to ranks $z_{(i)j}$ for each sample $j$ and normalize further $r_{ij}=|p/2-z_{(i)j}|$ to make the ranks symmetric around zero (Figure~\ref{methods}, step 2). This is done to up-weight the two tails of the rank distribution when computing the final enrichment score. We assess the enrichment score similar to the GSEA and ASSESS methods \cite{subramanian_gene_2005,edelman_analysis_2006} using the Kolmogorov-Smirnov (KS) like random walk statistic (Figure~\ref{methods}, step 3): \begin{equation} \label{walk} \nu_{jk}(\ell) = \frac{\sum_{i=1}^\ell |r_{ij}|^{\tau} I(g_{(i)} \in \gamma_k)}{\sum_{i=1}^p |r_{ij}|^{\tau}I(g_{(i)} \in \gamma_k)} - \frac{\sum_{i=1}^\ell I(g_{(i)} \not\in \gamma_k)}{p-|\gamma_k|}, \end{equation} where $\tau$ is a parameter describing the weight of the tail in the random walk (default $\tau = 1$), $\gamma_k$ is the $k$-th gene set, $I(g_{(i)} \in \gamma_k)$ is the indicator function on whether the $i$-th gene (the gene corresponding to the $i$-th ranked expression-level statistic) belongs to gene set $\gamma_k$, $|\gamma_k|$ is the number of genes in the $k$-th gene set, and $p$ is the number of genes in the data set. Conceptually, Eq.~\ref{walk} produces a distribution over the genes to assess if the genes in the gene set are more likely to be found at either tail of the rank distribution (see \cite{subramanian_gene_2005,edelman_analysis_2006} for a more detailed description). 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 \cite{subramanian_gene_2005,edelman_analysis_2006,verhaak_integrated_2010} and a normalized ES. The first ES is the maximum deviation from zero of the random walk of the $j$-th sample with respect to the $k$-th gene set : \begin{equation} \label{escore} ES^{\mbox{\tiny{max}}}_{jk} = \nu_{jk}[\arg \max_{\ell=1,\dots,p} \left|\nu_{jk}(\ell)\right|]. \end{equation} For each gene set $k$, this approach produces a distribution of enrichment scores that is bimodal (Figure~\ref{methods}, step 4, top panel). This is an intrinsic property of the KS like random walk, which generates non-zero maximum deviations under the null distribution. In GSEA \cite{subramanian_gene_2005} it is also observed that the empirical null distribution obtained by permuting phenotypes is bimodal and, for this reason, significance is determined independently using the positive and negative sides of the null distribution. In our case, we would like to provide a standard Gaussian distribution of enrichment scores under the null hypothesis of no change in pathway activity throughout the sample population. For this purpose we propose a second, alternative score that produces an ES distribution approximating this requirement (Figure~\ref{methods}, step 4, bottom panel): \begin{equation} \label{escore_2} ES^{\mbox{\tiny{diff}}}_{jk} = \left|ES^{+}_{jk}\right| - \left|ES^{-}_{jk}\right|=\max_{\ell=1,\dots,p}(0,\nu_{jk}(\ell)) - \min_{\ell=1,\dots,p}(0,\nu_{jk}(\ell))\,, \end{equation} where $ES_{jk}^{+}$ and $ES_{jk}^{-}$ are the largest positive and negative random walk deviations from zero, respectively, for sample $j$ and gene set $k$. This statistic may be compared to the Kuiper test statistic \cite{pearson_comparison_1963}, which sums the maximum and minimum deviations to make the test statistic more sensitive in the tails. In contrast, our test statistic penalizes deviations that are large in both tails, and provides a ``normalization'' of the enrichment score by subtracting potential noise. There is a clear biological interpretation of this statistic, it emphasizes genes in pathways that are concordantly activated in one direction only, either over-expressed or under-expressed relative to the overall population. For pathways containing genes strongly acting in both directions, the deviations will cancel each other out and show little or no enrichment. Because this statistic is unimodal and approximately normal, downstream analyses which may impose distributional assumptions on the data are possible. Figure~\ref{methods}, step 4 shows a simple simulation where standard Gaussian deviates are independenty sampled from $p=20,000$ genes and $n=30$ samples, thus mimicking a null distribution of no change in gene expression. One hundred gene sets are uniformly sampled at random from the $p$ genes with sizes ranging from 10 to 100 genes. Using these two inputs, we calculate the maximum deviation ES and the normalized ES. The resulting distributions are depicted in Figure~\ref{methods}, step 4 and in the larger figure below, illustrating the previous description. <<>>= 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.obs es.dif <- gsva(X, gs, mx.diff=TRUE, verbose=FALSE, parallel.sz=1)$es.obs @ \begin{center} <>= 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) @ \end{center} \bigskip Although the GSVA algorithm itself does not evaluate statistical significance for the enrichment of gene sets, significance with respect to one or more phenotypes can be easily evaluated using conventional statistical models. Likewise, false discovery rates can be estimated by permuting the sample labels (Methods). Examples of these techniques are provided in the following section. \section{Overview of the package} The \Rpackage{GSVA} package implements the methodology described in the previous section in the function \Rfunction{gsva()} which requires two main input arguments: the gene expression data and a collection of gene sets. The expression data can be provided either as a \Rclass{matrix} object of genes (rows) by sample (columns) expression values, or as an \Rclass{ExpressionSet} object. The collection of gene sets can be provided either as a \Rclass{list} object with names identifying gene sets and each entry of the list containing the gene identifiers of the genes forming the corresponding set, or as a \Rclass{GeneSetCollection} object as defined in the \Rpackage{GSEABase} package. When the two main arguments are an \Rclass{ExpressionSet} object and a \Rclass{GeneSetCollection} object, the \Rfunction{gsva()} function will first translate the gene identifiers used in the \Rclass{GeneSetCollection} object into the corresponding feature identifiers of the \Rclass{ExpressionSet} object, according to its corresponding annotation package. This translation is carried out by an internal call to the function \Rfunction{mapIdentifiers()} from the \Rpackage{GSEABase} package. This means that both input arguments may specify features with different types of identifiers, such as Entrez IDs and probeset IDs, and the \Rpackage{GSEABase} package will take care of mapping them to each other. A second filtering step is applied that removes genes without matching features in the \Rclass{ExpressionSet} object. If the expression data is given as a \Rclass{matrix} object then only the latter filtering step will be applied by the \Rfunction{gsva()} function and, therefore, it will be the responsibility of the user to have the same type of identifiers in both the expression data and the gene sets. After these automatic filtering steps, we may additionally filter out gene sets that do not meet a minimum and/or maximum size specified by the optional arguments \Robject{min.sz} and \Robject{max.sz} which are set by default to 1 and \Robject{Inf}, respectively. Finally, the \Rfunction{gsva()} function will carry out the calculations specified in the previous section and return a gene set by sample matrix of GSVA enrichment scores in the form of a \Rclass{matrix} object when this was the class of the input expression data object. Otherwise, it will return an \Rclass{ExpressionSet} object inheriting all the corresponding phenotypic information from the input data. An important argument of the \Rfunction{gsva()} function is the flag \Robject{mx.diff} which is set to \Robject{TRUE} by default. Under this default setting, GSVA enrichment scores are calculated using Equation~\ref{escore_2}, and therefore, are more amenable by analysis techniques that assume the data to be normally distributed. When setting \Robject{mx.diff=FALSE}, then Equation~\ref{escore} is employed, calculating enrichment in an analogous way to classical GSEA which typically provides a bimodal distribution of GSVA enrichment scores for each gene. Since the calculations for each gene set are independent from each other, the \Rfunction{gsva()} function offers two possibilities to perform them in parallel. One consists of loading the library \Rpackage{snow}, which will enable the parallelization of the calculations through a cluster of computers. In order to activate this option we should specify in the argument \Robject{parallel.sz} the number of processors we want to use (default is zero which means no parallelization is going to be employed). The other is loading the library \Rpackage{parallel} and then the \Rfunction{gsva()} function will use the core processors of the computer where R is running. If we want to limit \Rfunction{gsva()} in the number of core processors that should be used, we can do it by specifying the number of cores in the \Robject{parallel.sz} argument. The other two functions of the \Rpackage{GSVA} package are \Rfunction{filterGeneSets()} and \Rfunction{computeGeneSetsOverlaps()} that serve to explicitly filter gene sets by size and by pair-wise overlap, respectively. Note that the size filter can also be applied within the \Rfunction{gsva()} function call. The \Rfunction{gsva()} function also offers the following three other unsupervised GSE methods that calculate single sample pathway summaries of expression and which can be selected through the \Robject{method} argument: \begin{itemize} \item \Robject{method="plage"} \citep{tomfohr_pathway_2005}. Pathway level analysis of gene expression (PLAGE) standardizes first expression profiles into z-scores over the samples and then calculates the singular value decomposition $Z_\gamma=UDV'$ on the z-scores of the genes in the gene set. The coefficients of the first right-singular vector (first column of $V$) are taken as the gene set summaries of expression over the samples. \item \Robject{method="zscore"} \citep{lee_inferring_2008}. The combined z-score method also, as PLAGE, standardizes first expression profiles into z-scores over the samples, but combines them together for each gene set at each individual sample as follows. Given a gene set $\gamma=\{1,\dots,k\}$ with z-scores $Z_1,\dots,Z_k$ for each gene, the combined z-score $Z_\gamma$ for the gene set $\gamma$ is defined as: \begin{equation} Z_\gamma = \frac{\sum_{i=1}^k Z_i}{\sqrt{k}}\,. \end{equation} \item \Robject{method="ssgsea"} \citep{barbie_systematic_2009}. Single sample GSEA (ssGSEA) calculates a gene set enrichment score per sample as the normalized difference in empirical cumulative distribution functions of gene expression ranks inside and outside the gene set. \end{itemize} By default \Robject{method="gsva"} and the \Rfunction{gsva()} function uses the GSVA algorithm. \section{Applications} In this section we illustrate the following applications of \Rpackage{GSVA}: \begin{itemize} \item Functional enrichment between two subtypes of leukemia. \item Identification of molecular signatures in distinct glioblastoma subtypes. \end{itemize} Throughout this vignette we will use the C2 collection of curated gene sets that form part of the Molecular Signatures Database (MSigDB) version 3.0. This particular collection of gene sets is provided as a \Rclass{GeneSetCollection} object called \Robject{c2BroadSets} in the accompanying experimental data package \Rpackage{GSVAdata}, which stores these and other data employed in this vignette. These data can be loaded as follows: <>= library(GSEABase) library(GSVAdata) data(c2BroadSets) c2BroadSets @ where we observe that \Robject{c2BroadSets} contains \Sexpr{length(c2BroadSets)} gene sets. We also need to load the following additional libraries: <>= library(Biobase) library(genefilter) library(limma) library(RColorBrewer) library(GSVA) @ As a final setup step for this vignette, we will employ the \Rfunction{cache()} function from the \Rpackage{Biobase} package in order to load some pre-computed results and speed up the building time of the vignette: <<>>= cacheDir <- system.file("extdata", package="GSVA") cachePrefix <- "cache4vignette_" @ In order to enforce re-calculating everything, either the call to the \Rfunction{cache()} function should be replaced by its first argument, or the following command should be written in the R console at this point: <>= file.remove(paste(cacheDir, list.files(cacheDir, pattern=cachePrefix), sep="/")) @ \subsection{Functional enrichment} In this section we illustrate how to identify functionally enriched gene sets between two phenotypes. As in most of the applications we start by calculating GSVA enrichment scores and afterwards, we will employ the linear modeling techniques implemented in the \Rpackage{limma} package to find the enriched gene sets. The data set we use in this section corresponds to the microarray data from \citep{armstrong_mll_2002} which consists of 37 different individuals with human acute leukemia, where 20 of them have conventional childhood acute lymphoblastic leukemia (ALL) and the other 17 are affected with the MLL (mixed-lineage leukemia gene) translocation. This leukemia data set is stored as an \verb+ExpressionSet+ object called \Robject{leukemia} in the \Rpackage{GSVAdata} package and details on how the data was pre-processed can be found in the corresponding help page. Enclosed with the RMA expression values we provide some metadata including the main phenotype corresponding to the leukemia sample subtype. <<>>= data(leukemia) leukemia_eset head(pData(leukemia_eset)) table(leukemia_eset$subtype) @ Let's examine the variability of the expression profiles across samples by plotting the cumulative distribution of IQR values as shown in Figure~\ref{figIQR}. About 50\% of the probesets show very limited variability across samples and, therefore, in the following non-specific filtering step we remove this fraction from further analysis. <>= png(filename="GSVA-figIQR.png", width=500, height=500, res=150) IQRs <- esApply(leukemia_eset, 1, IQR) plot.ecdf(IQRs, pch=".", xlab="Interquartile range (IQR)", main="Leukemia data") abline(v=quantile(IQRs, prob=0.5), lwd=2, col="red") dev.off() @ \begin{figure}[ht] \centerline{\includegraphics[width=0.5\textwidth]{GSVA-figIQR}} \caption{Empirical cumulative distribution of the interquartile range (IQR) of expression values in the leukemia data. The vertical red bar is located at the 50\% quantile value of the cumulative distribution.} \label{figIQR} \end{figure} We carry out a non-specific filtering step by discarding the 50\% of the probesets with smaller variability, probesets without Entrez ID annotation, probesets whose associated Entrez ID is duplicated in the annotation, and Affymetrix quality control probes: <<>>= 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 @ The calculation of GSVA enrichment scores is performed in one single call to the \verb+gsva()+ function. However, one should take into account that this function performs further non-specific filtering steps prior to the actual calculations. On the one hand, it matches gene identifiers between gene sets and gene expression values. On the other hand, it discards gene sets that do not meet minimum and maximum gene set size requirements specified with the arguments \verb+min.sz+ and \verb+max.sz+, respectively, which, in the call below, are set to 10 and 500 genes. Because we want to use \Rpackage{limma} on the resulting GSVA enrichment scores, we leave deliberately unchanged the default argument \Robject{mx.diff=TRUE} to obtain approximately normally distributed ES. <<>>= cache(leukemia_es <- gsva(leukemia_filtered_eset, c2BroadSets, min.sz=10, max.sz=500, verbose=TRUE)$es.obs, dir=cacheDir, prefix=cachePrefix) @ We test whether there is a difference between the GSVA enrichment scores from each pair of phenotypes using a simple linear model and moderated t-statistics computed by the \verb+limma+ package using an empirical Bayes shrinkage method \citep[see][]{Smyth_2004}. We are going to examine both, changes at gene level and changes at pathway level and since, as we shall see below, there are plenty of them, we are going to employ the following stringent cut-offs to attain a high level of statistical and biological significance: <<>>= adjPvalueCutoff <- 0.001 logFCcutoff <- log2(2) @ where we will use the latter only for the gene-level differential expression analysis. <<>>= 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 \Sexpr{sum(res[, "MLLvsALL"]!=0)} MSigDB C2 curated pathways that are differentially activated between MLL and ALL at \Sexpr{adjPvalueCutoff*100}\% FDR. When we carry out the corresponding differential expression analysis at gene level: <<>>= 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, \Sexpr{sum(res[, "MLLvsALL"]!=0)} genes show up as being differentially expressed with a minimum fold-change of \Sexpr{2^logFCcutoff} at \Sexpr{adjPvalueCutoff*100}\% FDR. We illustrate the genes and pathways that are changing by means of volcano plots (Fig.~\ref{leukemiaVolcano}. <>= png(filename="GSVA-leukemiaVolcano.png", width=800, height=500) par(mfrow=c(1,2)) plot(allGeneSets$logFC, -log10(allGeneSets$P.Value), pch=".", cex=4, col=grey(0.75), main="Gene sets", xlab="GSVA enrichment score difference", ylab=expression(-log[10]~~~Raw~P-value)) abline(h=-log10(max(allGeneSets$P.Value[allGeneSets$adj.P.Val <= adjPvalueCutoff])), col=grey(0.5), lwd=1, lty=2) points(allGeneSets$logFC[match(rownames(DEgeneSets), rownames(allGeneSets))], -log10(allGeneSets$P.Value[match(rownames(DEgeneSets), rownames(allGeneSets))]), pch=".", cex=4, col="red") text(max(allGeneSets$logFC)*0.85, -log10(max(allGeneSets$P.Value[allGeneSets$adj.P.Val <= adjPvalueCutoff])), sprintf("%.1f%% FDR", 100*adjPvalueCutoff), pos=1) plot(allGenes$logFC, -log10(allGenes$P.Value), pch=".", cex=4, col=grey(0.75), main="Genes", xlab="Log fold-change", ylab=expression(-log[10]~~~Raw~P-value)) abline(h=-log10(max(allGenes$P.Value[allGenes$adj.P.Val <= adjPvalueCutoff])), col=grey(0.5), lwd=1, lty=2) abline(v=c(-logFCcutoff, logFCcutoff), col=grey(0.5), lwd=1, lty=2) points(allGenes$logFC[match(rownames(DEgenes), rownames(allGenes))], -log10(allGenes$P.Value[match(rownames(DEgenes), rownames(allGenes))]), pch=".", cex=4, col="red") text(max(allGenes$logFC)*0.85, -log10(max(allGenes$P.Value[allGenes$adj.P.Val <= adjPvalueCutoff])), sprintf("%.1f%% FDR", 100*adjPvalueCutoff), pos=1) dev.off() @ \begin{figure} \centerline{\includegraphics[width=0.8\textwidth]{GSVA-leukemiaVolcano}} \caption{Volcano plots for differential pathway activation (left) and differential gene expression (right) in the leukemia data set.} \label{leukemiaVolcano} \end{figure} The signatures of both, the differentially activated pathways reported by the GSVA analysis and of the differentially expressed genes are shown in Figures \ref{leukemiaHeatmapGeneSets} and \ref{leukemiaHeatmapGenes}, respectively. Many of the gene sets and pathways reported in Figure~\ref{leukemiaHeatmapGeneSets} are directly related to ALL and MLL. <>= png(filename="GSVA-leukemiaHeatmapGeneSets.png", width=500, height=500) GSVAsco <- exprs(leukemia_es[rownames(DEgeneSets), ]) colorLegend <- c("darkred", "darkblue") names(colorLegend) <- c("ALL", "MLL") sample.color.map <- colorLegend[pData(leukemia_es)[, "subtype"]] names(sample.color.map) <- colnames(GSVAsco) sampleClustering <- hclust(as.dist(1-cor(GSVAsco, method="spearman")), method="complete") geneSetClustering <- hclust(as.dist(1-cor(t(GSVAsco), method="pearson")), method="complete") heatmap(GSVAsco, ColSideColors=sample.color.map, xlab="samples", ylab="Gene sets and pathways", margins=c(2, 20), labRow=substr(gsub("_", " ", gsub("^KEGG_|^REACTOME_|^BIOCARTA_", "", rownames(GSVAsco))), 1, 35), labCol="", scale="row", Colv=as.dendrogram(sampleClustering), Rowv=as.dendrogram(geneSetClustering)) legend("topleft", names(colorLegend), fill=colorLegend, inset=0.01, bg="white") dev.off() @ \begin{figure}[ht] \centerline{\includegraphics[width=0.7\textwidth]{GSVA-leukemiaHeatmapGeneSets}} \caption{Heatmap of differentially activated pathways at \Sexpr{adjPvalueCutoff*100}\% FDR in the Leukemia data set.} \label{leukemiaHeatmapGeneSets} \end{figure} <>= png(filename="GSVA-leukemiaHeatmapGenes.png", width=500, height=500) exps <- exprs(leukemia_eset[rownames(DEgenes), ]) colorLegend <- c("darkred", "darkblue") names(colorLegend) <- c("ALL", "MLL") sample.color.map <- colorLegend[pData(leukemia_eset)[, "subtype"]] names(sample.color.map) <- colnames(exps) sampleClustering <- hclust(as.dist(1-cor(exps, method="spearman")), method="complete") geneClustering <- hclust(as.dist(1-cor(t(exps), method="pearson")), method="complete") heatmap(exps, ColSideColors=sample.color.map, xlab="samples", ylab="Genes", labRow="", labCol="", scale="row", Colv=as.dendrogram(sampleClustering), Rowv=as.dendrogram(geneClustering), margins=c(2,2)) legend("topleft", names(colorLegend), fill=colorLegend, inset=0.01, bg="white") dev.off() @ \begin{figure}[ht] \centerline{\includegraphics[width=0.5\textwidth]{GSVA-leukemiaHeatmapGenes}} \caption{Heatmap of differentially expressed genes with a minimum fold-change of \Sexpr{2^logFCcutoff} at \Sexpr{adjPvalueCutoff*100}\% FDR in the leukemia data set.} \label{leukemiaHeatmapGenes} \end{figure} \subsection{Molecular signature identification} In \citep{verhaak_integrated_2010} four subtypes of Glioblastoma multiforme (GBM) - proneural, classical, neural and mesenchymal - were identified by the characterization of distinct gene-level expression patterns. Using eight gene set signatures specific to brain cell types - astrocytes, oligodendrocytes, neurons and cultured astroglial cells - derived from murine models by \citep{cahoy_transcriptome_2008}, we replicate the analysis of \citep{verhaak_integrated_2010} by employing GSVA to transform the gene expression measurements into enrichment scores for these eight gene sets, without taking the sample subtype grouping into account. We start by loading and have a quick glance to the data which forms part of the \verb+GSVAdata+ package: <<>>= data(gbm_VerhaakEtAl) gbm_eset head(featureNames(gbm_eset)) table(gbm_eset$subtype) data(brainTxDbSets) sapply(brainTxDbSets, length) lapply(brainTxDbSets, head) @ GSVA enrichment scores for the gene sets contained in \Robject{brainTxDbSets} are calculated, in this case using \Robject{mx.diff=FALSE}, as follows: <<>>= gbm_es <- gsva(gbm_eset, brainTxDbSets, mx.diff=FALSE, verbose=FALSE, parallel.sz=1)$es.obs @ Figure \ref{gbmSignature} shows the GSVA enrichment scores obtained for the up-regulated gene sets across the samples of the four GBM subtypes. As expected, the \emph{neural} class is associated with the neural gene set and the astrocytic gene sets. The \emph{mesenchymal} subtype is characterized by the expression of mesenchymal and microglial markers, thus we expect it to correlate with the astroglial gene set. The \emph{proneural} subtype shows high expression of oligodendrocytic development genes, thus it is not surprising that the oligodendrocytic gene set is highly enriched for ths group. Interestingly, the \emph{classical} group correlates highly with the astrocytic gene set. In summary, the resulting GSVA enrichment scores recapitulate accurately the molecular signatures from \citet{verhaak_integrated_2010}. <>= png(filename="GSVA-gbmSignature.png", width=700, height=500) subtypeOrder <- c("Proneural", "Neural", "Classical", "Mesenchymal") sampleOrderBySubtype <- sort(match(gbm_es$subtype, subtypeOrder), index.return=TRUE)$ix subtypeXtable <- table(gbm_es$subtype) subtypeColorLegend <- c(Proneural="red", Neural="green", Classical="blue", Mesenchymal="orange") geneSetOrder <- c("astroglia_up", "astrocytic_up", "neuronal_up", "oligodendrocytic_up") geneSetLabels <- gsub("_", " ", geneSetOrder) hmcol <- colorRampPalette(brewer.pal(10, "RdBu"))(256) hmcol <- hmcol[length(hmcol):1] heatmap(exprs(gbm_es)[geneSetOrder, sampleOrderBySubtype], Rowv=NA, Colv=NA, scale="row", margins=c(3,5), col=hmcol, ColSideColors=rep(subtypeColorLegend[subtypeOrder], times=subtypeXtable[subtypeOrder]), labCol="", gbm_es$subtype[sampleOrderBySubtype], labRow=paste(toupper(substring(geneSetLabels, 1,1)), substring(geneSetLabels, 2), sep=""), cexRow=2, main=" \n ") par(xpd=TRUE) text(0.22,1.11, "Proneural", col="red", cex=1.2) text(0.36,1.11, "Neural", col="green", cex=1.2) text(0.48,1.11, "Classical", col="blue", cex=1.2) text(0.66,1.11, "Mesenchymal", col="orange", cex=1.2) mtext("Gene sets", side=4, line=0, cex=1.5) mtext("Samples ", side=1, line=4, cex=1.5) dev.off() @ \begin{figure} \centerline{\includegraphics[width=0.6\textwidth]{GSVA-gbmSignature}} \caption{Heatmap of GSVA scores for cell-type brain signatures from murine models (y-axis) across GBM samples grouped by GBM subtype.} \label{gbmSignature} \end{figure} \section{Comparison with other methods} In this section we compare with simulated data the performance of GSVA with other methods producing pathway summaries of gene expression, concretely, PLAGE, the combined z-score and ssGSEA which are available through the argument \Robject{method} of the function \Rfunction{gsva()}. We employ the following simple linear additive model for simulating normalized microarray data on $p$ genes and $n$ samples divided in two groups representing a case-control scenario: \begin{equation} y_{ij} = \alpha_i + \beta_j + \epsilon_{ij}\,, \end{equation} where $\alpha_i\sim\mathcal{N}(0, 1)$ is a gene-specific effect, such as a probe-effect, with $i=1,\dots,p$, $\beta_j\sim\mathcal{N}(\mu_j, 0.5)$ is a sample-effect with $j=1,2$ and $e_{ij}\sim\mathcal{N}(0,1)$ corresponds to random noise. We will assess the statistical power to detect one differentially expressed (DE) gene set formed by 30 genes, out of $p=1000$, as a function of the sample size and two varying conditions: the fraction of differentially expressed genes in the gene set (50\% and 80\%) and the signal-to-noise ratio expressed as the magnitude of the mean sample effect for one of the sample groups ($\mu_1=0$ and either $\mu_2=0.5$ or $\mu_2=1$). Simulatenously, we will assess the empirical type-I error rate by building using one gene set of 30 non-DE genes. The following function enables simulating such data, computes the corresponding GSE scores with each method, performs a $t$-test on the tested gene sets between the two groups of samples for each method and returns the corresponding $p$-values: <<>>= runSim <- function(p, n, gs.sz, S2N, fracDEgs) { sizeDEgs <- round(fracDEgs * gs.sz) group.n <- round(n / 2) sampleEffect <- rnorm(n, mean=0, sd=1) sampleEffectDE <- rnorm(n, mean=S2N, sd=0.5) probeEffect <- rnorm(p, mean=0, sd=1) noise <- matrix(rnorm(p*n, mean=0, sd=1), nrow=p, ncol=n) noiseDE <- matrix(rnorm(p*n, mean=0, sd=1), nrow=p, ncol=n) M <- outer(probeEffect, sampleEffect, "+") + noise M2 <- outer(probeEffect, sampleEffectDE, "+") + noiseDE M[1:sizeDEgs, 1:group.n] <- M2[1:sizeDEgs, 1:group.n] rownames(M) <- paste0("g", 1:nrow(M)) geneSets <- list(H1GeneSet=paste0("g", 1:(gs.sz)), H0GeneSet=paste0("g", (gs.sz+1):(2*gs.sz))) es.gsva <- gsva(M, geneSets, verbose=FALSE, parallel.sz=1)$es.obs es.ss <- gsva(M, geneSets, method="ssgsea", verbose=FALSE, parallel.sz=1) es.z <- gsva(M, geneSets, method="zscore", verbose=FALSE, parallel.sz=1) es.plage <- gsva(M, geneSets, method="plage", verbose=FALSE, parallel.sz=1) h1.gsva.pval <- t.test(es.gsva["H1GeneSet", 1:group.n],es.gsva["H1GeneSet", (group.n+1):n])$p.value h1.ssgsea.pval <- t.test(es.ss["H1GeneSet", 1:group.n],es.ss["H1GeneSet", (group.n+1):n])$p.value h1.zscore.pval <- t.test(es.z["H1GeneSet", 1:group.n],es.z["H1GeneSet", (group.n+1):n])$p.value h1.plage.pval <- t.test(es.plage["H1GeneSet", 1:group.n],es.plage["H1GeneSet", (group.n+1):n])$p.value h0.gsva.pval <- t.test(es.gsva["H0GeneSet", 1:group.n],es.gsva["H0GeneSet", (group.n+1):n])$p.value h0.ssgsea.pval <- t.test(es.ss["H0GeneSet", 1:group.n],es.ss["H0GeneSet", (group.n+1):n])$p.value h0.zscore.pval <- t.test(es.z["H0GeneSet", 1:group.n],es.z["H0GeneSet", (group.n+1):n])$p.value h0.plage.pval <- t.test(es.plage["H0GeneSet", 1:group.n],es.plage["H0GeneSet", (group.n+1):n])$p.value c(h1.gsva.pval, h1.ssgsea.pval, h1.zscore.pval, h1.plage.pval, h0.gsva.pval, h0.ssgsea.pval, h0.zscore.pval, h0.plage.pval) } @ The next function takes the $p$-values of the output of the previous function and estimates the statistical power as the fraction of non-rejections on the DE gene set, and the empirical type-I error rate as the fraction of rejections on the non-DE gene set, at a significant level $\alpha=0.05$. <<>>= estPwrTypIerr <- function(pvals, alpha=0.05) { N <- ncol(pvals) c(1 - sum(pvals[1, ] > alpha)/N, 1 - sum(pvals[2, ] > alpha)/N,1 - sum(pvals[3, ] > alpha)/N, 1 - sum(pvals[4, ] > alpha)/N, sum(pvals[5, ] <= alpha)/N, sum(pvals[6, ] <= alpha)/N, sum(pvals[7, ] <= alpha)/N, sum(pvals[8, ] <= alpha)/N) } @ Finally, we perform the simulation on each of the four described scenarios 60 times using the code below. The results in Fig.~\ref{simpower} show that GSVA attains higher statistical power than the other three methods in each of the simulated scenarios, while providing a similar control of the type-I error rate. Note that the fluctuations of this latter estimate are due to the limited number of times we repeat the simulation; see \citet{haenzelmann_castelo_guinney_2013} for more stable estimates obtained with a much larger number of repeated simulations. <<>>= set.seed(1234) exp1 <- cbind(estPwrTypIerr(replicate(60, runSim(1000, 10, gs.sz=30, S2N=0.5, fracDEgs=0.5))), estPwrTypIerr(replicate(60, runSim(1000, 20, gs.sz=30, S2N=0.5, fracDEgs=0.5))), estPwrTypIerr(replicate(60, runSim(1000, 40, gs.sz=30, S2N=0.5, fracDEgs=0.5))), estPwrTypIerr(replicate(60, runSim(1000, 60, gs.sz=30, S2N=0.5, fracDEgs=0.5)))) exp2 <- cbind(estPwrTypIerr(replicate(60, runSim(1000, 10, gs.sz=30, S2N=1.0, fracDEgs=0.5))), estPwrTypIerr(replicate(60, runSim(1000, 20, gs.sz=30, S2N=1.0, fracDEgs=0.5))), estPwrTypIerr(replicate(60, runSim(1000, 40, gs.sz=30, S2N=1.0, fracDEgs=0.5))), estPwrTypIerr(replicate(60, runSim(1000, 60, gs.sz=30, S2N=1.0, fracDEgs=0.5)))) exp3 <- cbind(estPwrTypIerr(replicate(60, runSim(1000, 10, gs.sz=30, S2N=0.5, fracDEgs=0.8))), estPwrTypIerr(replicate(60, runSim(1000, 20, gs.sz=30, S2N=0.5, fracDEgs=0.8))), estPwrTypIerr(replicate(60, runSim(1000, 40, gs.sz=30, S2N=0.5, fracDEgs=0.8))), estPwrTypIerr(replicate(60, runSim(1000, 60, gs.sz=30, S2N=0.5, fracDEgs=0.8)))) exp4 <- cbind(estPwrTypIerr(replicate(60, runSim(1000, 10, gs.sz=30, S2N=1.0, fracDEgs=0.8))), estPwrTypIerr(replicate(60, runSim(1000, 20, gs.sz=30, S2N=1.0, fracDEgs=0.8))), estPwrTypIerr(replicate(60, runSim(1000, 40, gs.sz=30, S2N=1.0, fracDEgs=0.8))), estPwrTypIerr(replicate(60, runSim(1000, 60, gs.sz=30, S2N=1.0, fracDEgs=0.8)))) @ <>= plotPower <- function(statpower, main, legendposition="bottomright", ...) { plot(statpower[1,], ylim=c(0, 1.0), type="b", lwd=2, pch=1, main=main, col="blue", ylab="Statistcal Power", xlab="Sample Size", xaxt="n") lines(statpower[2,], col="red", type="b", lwd=2, pch=2) lines(statpower[3,], col="darkgreen", type="b", lwd=2, pch=3) lines(statpower[4,], col="lightgreen", type="b", lwd=2, pch=4) if (!is.null(legendposition)) legend(legendposition, c("GSVA","ssGSEA","z-score","PLAGE"), col=c("blue","red","darkgreen","lightgreen"),pch=1:4,lty=1,lwd=2,inset=0.02) axis(1,at=1:4, labels=c("10","20","40","60")) } plotType1Error <- function(tmp, title, legendposition="bottomright", alpha=0.05, ...){ plot(tmp[5,],ylim=c(0, 0.2),type="b",lwd=2,pch=1, col="blue",ylab="Empirical Type-I Error",xlab="Sample Size",xaxt="n",main=title, ...) lines(tmp[6,],col="red",type="b",lwd=2,pch=2) lines(tmp[7,],col="darkgreen",type="b",lwd=2,pch=3) lines(tmp[8,],col="lightgreen",type="b",lwd=2,pch=4) if (!is.null(legendposition)) legend(legendposition,c("GSVA","ssGSEA","z-score","PLAGE"),col=c("blue","red","darkgreen","lightgreen"),pch=1:4,lty=1,lwd=2,inset=0.02) axis(1,at=c(1:dim(tmp)[2]), labels=c("10","20","40","60")) abline(h=alpha, lty=2) } labelPlot <- function(lab, font, cex, offsetx=0.05, offsety=0.05) { par(xpd=TRUE) w <- par("usr")[2] - par("usr")[1] h <- par("usr")[4] - par("usr")[3] text(par("usr")[1]-w*offsetx, par("usr")[4]+h*offsety, lab, font=font, cex=cex) par(xpd=FALSE) } par(mfrow=c(4,2), mar=c(4, 4, 2, 1)) plotPower(exp1, main="", legendposition=NULL, las=1) labelPlot("A", 2, 2, 0.2, 0.15) plotType1Error(exp1,"",legendposition="topright", las=1) labelPlot("B", 2, 2, 0.2, 0.15) plotPower(exp2, main="", legendposition=NULL, las=1) labelPlot("C", 2, 2, 0.2, 0.15) plotType1Error(exp2,"",legendposition="topright", las=1) labelPlot("D", 2, 2, 0.2, 0.15) plotPower(exp3, main="", legendposition=NULL, las=1) labelPlot("E", 2, 2, 0.2, 0.15) plotType1Error(exp3,"",legendposition="topright", las=1) labelPlot("F", 2, 2, 0.2, 0.15) plotPower(exp4, main="", legendposition=NULL, las=1) labelPlot("G", 2, 2, 0.2, 0.15) plotType1Error(exp4,"",legendposition="topright", las=1) labelPlot("H", 2, 2, 0.2, 0.15) @ \begin{figure}[p] \centerline{\includegraphics[height=0.8\textheight]{GSVA-powertype1errsim}} \caption{{\bf Comparison of the statistical power and type-I error rate between GSVA, PLAGE, single sample GSEA (ssGSEA) and combined z-score (zscore).} The averaged results of 60 simulations are depicted as function of the sample size on the $x$-axis, for each of the GSE methods. On the $y$-axis either the statistical power (A, C, E, G) or the empirical type-I error rate (B, D, F, H) is shown. GSE scores were calculated with each method with respect to two gene sets, one of them differentially expressed (DE) and the other one not. Statistical power and empirical type-I error rates were estimated by performing a $t$-test on the DE and non-DE gene sets, respectively, at a significance level of $\alpha=0.05$. These simulations were carried out under the following four different scenarios for the DE gene set: (A,B) weak signal-to-noise ratio, 50\% of DE genes in the DE gene set; (C,D) strong signal-to-noise ratio, 50\% of DE genes in the DE gene set; (E, F) weak signal-to-noise ratio, 80\% of DE genes in the DE gene set; (G, H) strong signal-to-noise ratio, 80\% of DE gene in the DE gene set.} \label{simpower} \end{figure} \section{GSVA for RNA-Seq data} In this section we illustrate how to use GSVA with RNA-seq data and, more importantly, how the method provides pathway activity profiles analogous to the ones obtained from microarray data by using samples of lymphoblastoid cell lines (LCL) from HapMap individuals which have been profiled using both technologies \cite{huang_genome-wide_2007, pickrell_understanding_2010}. These data form part of the experimental package \Rpackage{GSVAdata} and the corresponding help pages contain details on how the data were processed. We start loading these data and verifying that they indeed contain expression data for the same genes and samples, as follows: <<>>= data(commonPickrellHuang) stopifnot(identical(featureNames(huangArrayRMAnoBatchCommon_eset), featureNames(pickrellCountsArgonneCQNcommon_eset))) stopifnot(identical(sampleNames(huangArrayRMAnoBatchCommon_eset), sampleNames(pickrellCountsArgonneCQNcommon_eset))) @ Next, for the current analysis we use the subset of canonical pathways from the C2 collection of MSigDB Gene Sets. These correspond to the following pathways from KEGG, REACTOME and BIOCARTA: <<>>= canonicalC2BroadSets <- c2BroadSets[c(grep("^KEGG", names(c2BroadSets)), grep("^REACTOME", names(c2BroadSets)), grep("^BIOCARTA", names(c2BroadSets)))] canonicalC2BroadSets @ Additionally, we extend this collection of gene sets with two formed by genes with sex-specific expression: <<<>>= 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 @ We calculate now GSVA enrichment scores for these gene sets using first the microarray data and then the RNA-seq data. Note that the only requirement to do the latter is to set the argument \Robject{rnaseq=TRUE} which is \Robject{FALSE} by default. <<<>>= esmicro <- gsva(huangArrayRMAnoBatchCommon_eset, canonicalC2BroadSets, min.sz=5, max.sz=500, mx.diff=TRUE, verbose=FALSE, rnaseq=FALSE, parallel.sz=1)$es.obs dim(esmicro) esrnaseq <- gsva(pickrellCountsArgonneCQNcommon_eset, canonicalC2BroadSets, min.sz=5, max.sz=500, mx.diff=TRUE, verbose=FALSE, rnaseq=TRUE, parallel.sz=1)$es.obs dim(esrnaseq) @ To compare expression values from both technologies we are going to transform the RNA-seq read counts into RPKM values. For this purpose we need gene length and G+C content information also stored in the \Rpackage{GSVAdata} package and use the \Rfunction{cpm()} function from the \Rpackage{edgeR} package. Note that RPKMs can only be calculated for those genes for which the gene length and G+C content information is available: <<>>= library(edgeR) data(annotEntrez220212) head(annotEntrez220212) cpm <- cpm(exprs(pickrellCountsArgonneCQNcommon_eset)) dim(cpm) common <- intersect(rownames(cpm), rownames(annotEntrez220212)) length(common) rpkm <- sweep(cpm[common, ], 1, annotEntrez220212[common, "Length"] / 10^3, FUN="/") dim(rpkm) dim(huangArrayRMAnoBatchCommon_eset[rownames(rpkm), ]) @ We finally calculate Spearman correlations between gene and gene-level expression values and gene set level GSVA enrichment scores produced from data obtained by microarray and RNA-seq technologies: <<>>= corsrowsgene <- sapply(1:nrow(huangArrayRMAnoBatchCommon_eset[rownames(rpkm), ]), function(i, expmicro, exprnaseq) cor(expmicro[i, ], exprnaseq[i, ], method="pearson"), exprs(huangArrayRMAnoBatchCommon_eset[rownames(rpkm), ]), log2(rpkm+0.1)) names(corsrowsgene) <- rownames(rpkm) corsrowsgs <- sapply(1:nrow(esmicro), function(i, esmicro, esrnaseq) cor(esmicro[i, ], esrnaseq[i, ], method="spearman"), exprs(esmicro), exprs(esrnaseq)) names(corsrowsgs) <- rownames(esmicro) @ In panels A and B of Figure~\ref{rnaseqcomp} we can see the distribution of these correlations at gene and gene set level. They show that GSVA enrichment scores correlate similarly to gene expression levels produced by both profiling technologies. <>= png(filename="GSVA-RNAseqComp.png", width=1100, height=1100, res=150) par(mfrow=c(2,2), mar=c(4, 5, 3, 2)) hist(corsrowsgene, xlab="Spearman correlation", main="Gene level\n(RNA-seq RPKM vs Microarray RMA)", xlim=c(-1, 1), col="grey", las=1) par(xpd=TRUE) text(par("usr")[1]*1.5, par("usr")[4]*1.1, "A", font=2, cex=2) par(xpd=FALSE) hist(corsrowsgs, xlab="Spearman correlation", main="Gene set level\n(GSVA enrichment scores)", xlim=c(-1, 1), col="grey", las=1) par(xpd=TRUE) text(par("usr")[1]*1.5, par("usr")[4]*1.1, "B", font=2, cex=2) par(xpd=FALSE) plot(exprs(esrnaseq)["MSY", ], exprs(esmicro)["MSY", ], xlab="GSVA scores RNA-seq", ylab="GSVA scores microarray", main=sprintf("MSY R=%.2f", cor(exprs(esrnaseq)["MSY", ], exprs(esmicro)["MSY", ])), las=1, type="n") sprintf("MSY R=%.2f", cor(exprs(esrnaseq)["MSY", ], exprs(esmicro)["MSY", ])) abline(lm(exprs(esmicro)["MSY", ] ~ exprs(esrnaseq)["MSY", ]), lwd=2, lty=2, col="grey") points(exprs(esrnaseq)["MSY", pickrellCountsArgonneCQNcommon_eset$Gender == "Female"], exprs(esmicro)["MSY", huangArrayRMAnoBatchCommon_eset$Gender == "Female"], col="red", pch=21, bg="red", cex=1) points(exprs(esrnaseq)["MSY", pickrellCountsArgonneCQNcommon_eset$Gender == "Male"], exprs(esmicro)["MSY", huangArrayRMAnoBatchCommon_eset$Gender == "Male"], col="blue", pch=21, bg="blue", cex=1) par(xpd=TRUE) text(par("usr")[1]*1.5, par("usr")[4]*1.1, "C", font=2, cex=2) par(xpd=FALSE) plot(exprs(esrnaseq)["XiE", ], exprs(esmicro)["XiE", ], xlab="GSVA scores RNA-seq", ylab="GSVA scores microarray", main=sprintf("XiE R=%.2f", cor(exprs(esrnaseq)["XiE", ], exprs(esmicro)["XiE", ])), las=1, type="n") abline(lm(exprs(esmicro)["XiE", ] ~ exprs(esrnaseq)["XiE", ]), lwd=2, lty=2, col="grey") points(exprs(esrnaseq["XiE", pickrellCountsArgonneCQNcommon_eset$Gender == "Female"]), exprs(esmicro)["XiE", huangArrayRMAnoBatchCommon_eset$Gender == "Female"], col="red", pch=21, bg="red", cex=1) points(exprs(esrnaseq)["XiE", pickrellCountsArgonneCQNcommon_eset$Gender == "Male"], exprs(esmicro)["XiE", huangArrayRMAnoBatchCommon_eset$Gender == "Male"], col="blue", pch=21, bg="blue", cex=1) par(xpd=TRUE) text(par("usr")[1]*1.5, par("usr")[4]*1.1, "D", font=2, cex=2) par(xpd=FALSE) dev.off() @ \begin{figure}[p] \centerline{\includegraphics[height=0.5\textheight]{GSVA-RNAseqComp}} \caption{{\bf GSVA for RNA-seq (Argonne).} A. Distribution of Spearman correlation values between gene expression profiles of RNA-seq and microarray data. B. Distribution of Spearman correlation values between GSVA enrichment scores of gene sets calculated from RNA-seq and microarray data. C and D. Comparison of GSVA enrichment scores obtained from microarray and RNA-seq data for two gene sets containing genes with sex-specific expression: MSY formed by genes from the male-specific region of the Y chromosome, thus male-specific, and XiE formed by genes that escape X-inactivation in females, thus female-specific. Red and blue dots represent female and male samples, respectively. In both cases GSVA-scores show very high correlation between the two profiling technologies where female samples show higher enrichment scores in the female-specific gene set and male samples show higher enrichment scores in the male-specific gene set.} \label{rnaseqcomp} \end{figure} We also examined the two gene sets containing gender specific genes in detail: those that escape X-inactivation in female samples \citep{carrel_x-inactivation_2005} and those that are located on the male-specific region of the Y chrosomome \citep{skaletsky_male-specific_2003}. In panels C and D of Figure~\ref{rnaseqcomp} we can see how microarray and RNA-seq enrichment scores correlate very well in these gene sets, with $\rho=0.82$ for the male-specific gene set and $\rho=0.78$ for the female-specific gene set. Male and female samples show higher GSVA enrichment scores in their corresponding gene set. This demonstrates the flexibility of GSVA to enable analogous unsupervised and single sample GSE analyses in data coming from both, microarray and RNA-seq technologies. \section{Session Information} <>= toLatex(sessionInfo()) @ \bibliography{GSVA} \bibliographystyle{apalike} \end{document}