@@ -1,24 +1,21 @@

 Package: GSVA

-Version: 1.39.15

+Version: 1.39.16

 Title: Gene Set Variation Analysis for microarray and RNA-seq data

 Authors@R: c(person("Justin", "Guinney", role=c("aut", "cre"), email="justin.guinney@sagebase.org"),

              person("Robert", "Castelo", role="aut", email="robert.castelo@upf.edu"),

              person("Alexey", "Sergushichev", role="ctb", email="alsergbox@gmail.com"),

              person("Pablo Sebastian", "Rodriguez", role="ctb", email="pablosebastian.rodriguez@upf.edu"))

 Depends: R (>= 3.5.0)

-Imports: methods, stats, utils, graphics, S4Vectors, IRanges,

-         Biobase, SummarizedExperiment, GSEABase, Matrix,

-         parallel, BiocParallel, SingleCellExperiment, 

-         sparseMatrixStats, DelayedArray, DelayedMatrixStats,

-         HDF5Array, BiocSingular

-Suggests: 

-    RUnit, BiocGenerics, limma, RColorBrewer, genefilter, edgeR, GSVAdata,

-    shiny, shinythemes, ggplot2, data.table, plotly

+Imports: methods, stats, utils, graphics, BiocGenerics, S4Vectors, IRanges,

+         Biobase, SummarizedExperiment, GSEABase, Matrix, parallel,

+         BiocParallel, SingleCellExperiment, sparseMatrixStats, DelayedArray,

+         DelayedMatrixStats, HDF5Array, BiocSingular

+Suggests: RUnit, BiocStyle, knitr, markdown, limma, RColorBrewer, genefilter,

+          edgeR, GSVAdata, shiny, shinythemes, ggplot2, data.table, plotly

 Description: Gene Set Variation Analysis (GSVA) is a non-parametric, unsupervised method for estimating variation of gene set enrichment through the samples of a expression data set. 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 the evaluation of pathway enrichment for each sample. This new matrix of GSVA enrichment scores facilitates applying standard analytical methods like functional enrichment, survival analysis, clustering, CNV-pathway analysis or cross-tissue pathway analysis, in a pathway-centric manner.

 License: GPL (>= 2)

-LazyLoad: yes

+VignetteBuilder: knitr

 URL: https://github.com/rcastelo/GSVA

 BugReports: https://github.com/rcastelo/GSVA/issues

 Encoding: latin1

-biocViews: Microarray, Pathways, GeneSetEnrichment

-Config/testthat/edition: 3

+biocViews: FunctionalGenomics, Microarray, RNASeq, Pathways, GeneSetEnrichment

@@ -978,14 +978,13 @@ setGeneric("filterGeneSets", function(gSets, ...) standardGeneric("filterGeneSet

 setMethod("filterGeneSets", signature(gSets="list"),

           function(gSets, min.sz=1, max.sz=Inf) {

-	gSetsLen <- sapply(gSets,length)

+	gSetsLen <- lengths(gSets)

 	return (gSets[gSetsLen >= min.sz & gSetsLen <= max.sz])	

 })

 setMethod("filterGeneSets", signature(gSets="GeneSetCollection"),

           function(gSets, min.sz=1, max.sz=Inf) {

-	gSetsLen <- sapply(geneIds(gSets),length)

-	return (gSets[gSetsLen >= min.sz & gSetsLen <= max.sz])	

+  filterGeneSets(geneIds(gSets), min.sz., max.sz)

 })

new file mode 100644

@@ -0,0 +1,340 @@

+---

+title: "GSVA: gene set variation analysis for bulk expression data"

+author:

+- name: Robert Castelo

+  affiliation:

+  - &idupf Dept. of Experimental and Health Sciences, Universitat Pompeu Fabra, Barcelona, Spain

+  email: robert.castelo@upf.edu

+- name: Pablo Sebastian Rodriguez

+  affiliation: *idupf

+  email: pablosebastian.rodriguez@upf.edu

+- name: Justin Guinney 

+  affiliation: 

+  - Sage Bionetworks

+  email: justin.guinney@sagebase.org

+abstract: >

+  The GSVA package provides the implementation of four single-sample gene set

+  enrichment methods, concretely _zscore_, _plage_, _ssGSEA_ and its own called

+  _GSVA_. These methods transform an input gene-by-sample expression data matrix

+  into a gene-set-by-sample expression data matrix. Thereby enabling the

+  estimation of pathway activity for each sample and facilitating pathway-centric

+  analyses of gene expression data. In this vignette we illustrate how to use

+  the GSVA package with bulk microarray and RNA-seq expression data.

+date: "`r BiocStyle::doc_date()`"

+package: "`r pkg_ver('GSVA')`"

+vignette: >

+  %\VignetteIndexEntry{GSVA for bulk expression data}

+  %\VignetteEngine{knitr::rmarkdown}

+  %\VignetteEncoding{UTF-8}

+  %\VignetteKeywords{GeneExpression, Microarray, RNAseq, GeneSetEnrichment, Pathway}

+output:

+  BiocStyle::html_document:

+    toc: true

+    toc_float: true

+    number_sections: true

+bibliography: GSVA.bib

+---

+

+**License**: `r packageDescription("GSVA")[["License"]]`

+

+```{r setup, include=FALSE}

+options(width=80)

+knitr::opts_chunk$set(collapse=TRUE,

+                      message=FALSE)

+```

+

+# Quick start

+

+`r Biocpkg("GSVA")` is an R package distributed as part of the

+[Bioconductor](https://bioconductor.org) project. To install the package, start

+R and enter:

+

+```{r library_install, message=FALSE, cache=FALSE, eval=FALSE}

+install.packages("BiocManager")

+BiocManager::install("GSVA")

+```

+

+Once `r Biocpkg("GSVA")` is installed, it can be loaded with the following command.

+

+```{r load_library, message=FALSE, warning=FALSE, cache=FALSE}

+library(GSVA)

+```

+

+Given a gene expression data matrix with rows corresponding to genes and columns

+to samples, such as this one simulated from random Gaussian data:

+

+```{r}

+p <- 10000 ## number of genes

+n <- 30    ## number of samples

+## simulate expression values from a standard Gaussian distribution

+X <- matrix(rnorm(p*n), nrow=p,

+            dimnames=list(paste0("g", 1:p), paste0("s", 1:n)))

+X[1:5, 1:5]

+```

+

+Given a collection of gene sets stored, for instance, in a `list` object such as

+this one with gene sampled uniformly at random without replacement into the gene sets:

+

+```{r}

+## sample gene set sizes

+gs <- as.list(sample(10:100, size=100, replace=TRUE))

+## sample gene sets

+gs <- lapply(gs, function(n, p)

+                   paste0("g", sample(1:p, size=n, replace=FALSE)), p)

+names(gs) <- paste0("gs", 1:length(gs))

+```

+

+We can calculate GSVA enrichment scores as follows:

+

+```{r}

+gsva.es <- gsva(X, gs, verbose=FALSE)

+dim(gsva.es)

+gsva.es[1:5, 1:5]

+```

+

+So, 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

+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 using GSVA on large molecular data sets and speed up computing time. You will

+find information on all these features in this vignette.

+

+# Introduction

+

+Gene set variation analysis (GSVA) provides an estimate of pathway activity

+by transforming an input gene-by-sample expression data matrix

+into a gene-set-by-sample one. 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).

+

+The GSVA package provides an implementation of this approach for the following

+methods:

+

+* _plage_ [@tomfohr_pathway_2005]. Pathway level analysis of gene expression

+  (PLAGE) standardizes expression profiles over the samples and then, for each

+  gene set, it performs a singular value decomposition (SVD) over its genes.

+  The coefficients of the first right-singular vector are returned as the

+  estimates of pathway activity over the samples. Note that, because of how

+  SVD is calculated, the sign of its singular vectors is arbitrary.

+

+* _zscore_ [@lee_inferring_2008]. The z-score method standardizes expression

+  profiles over the samples and then, for each gene set, combines the

+  standardized values as follows. Given a gene set $\gamma=\{1,\dots,k\}$

+  with standardized values $z_1,\dots,z_k$ for each gene in a specific sample,

+  the combined z-score $Z_\gamma$ for the gene set $\gamma$ is defined as:

+  $$

+  Z_\gamma = \frac{\sum_{i=1}^k z_i}{\sqrt{k}}\,.

+  $$

+

+* _ssgsea_ [@barbie_systematic_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_systematic_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.

+

+* _gsva_ [@haenzelmann_castelo_guinney_2013]. This is the default method of

+  the package and similarly to ssGSEA, is a non-parametric method that

+  uses the empirical CDFs of gene expression ranks inside and outside the gene

+  set, but it starts by calculating an expression-level statistic that brings

+  gene expression profiles with a different dynamic range to a common scale.

+

+The interested user may find full technical details about how these methods

+work in their corresponding article cited above. If you use any of them in a

+publication, please cite it with the given bibliographic reference.

+

+# Overview of the GSVA functionality

+

+The workhorse of the GSVA package is the function `gsva()`, which requires

+the following two input arguments:

+

+1. A normalized 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 `r Biocpkg("Biobase")`.

+   * A `SummarizedExperiment` object, see package

+     `r Biocpkg("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 `r Biocpkg("GSEABase")`.

+

+One advantage of providing the input data using specialized containers such as

+`ExpressionSet`, `SummarizedExperiment` and `GeneSetCollection` is that the

+`gsva()` function will automatically map the gene identifiers between the

+expression data and the gene sets (internally calling the function

+`mapIdentifiers()` from the package `r Biocpkg("GSEABase")`), when they come

+from different standard nomenclatures, i.e., _Ensembl_ versus _Entrez_, provided

+the input objects contain the appropriate metadata; see next section.

+

+If either the input gene expression data is provided as a `matrix` object or

+the gene sets are provided in a `list` object, or both, it is then the

+responsibility of the user to ensure that both objects contain gene identifiers

+following the same standard nomenclature.

+

+Before the actual calculations take place, the `gsva()` function will apply

+the following filters:

+

+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 filter, do not meet a

+   minimum and maximum size, which by default is one for the minimum size and

+   has no limit in the maximum size.

+

+If as a result of this filter either no genes or gene sets are left, the

+`gsva()` function will prompt an error. A common cause for an error at this

+stage is that gene identifiers between the expression data matrix and the gene

+sets do not belong to the same standard nomenclature and could not be mapped.

+

+By default, the `gsva()` function employs the method described by

+@haenzelmann_castelo_guinney_2013 but this can be changed using the argument

+`method`, which can take values `gsva` (default), `zscore`, `plage` or

+`ssgsea`, corresponding to the methods briefly described in the introduction.

+

+When `method="gsva"` (default), the user can additionally tune the following

+parameters:

+

+* `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 whih the largest positive and

+  negative random walk devations add 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`.

+

+In general, the default values for the previous parameters are suitable for

+most analysis settings, which usually consist of normalized continuous

+expression values.

+

+# Gene sets definitions and mapping to gene identifiers

+

+# Example applications

+

+## Molecular signature identification

+

+In [@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 four

+gene set signatures specific to brain cell types (astrocytes, oligodendrocytes,

+neurons and cultured astroglial cells), derived from murine models by

+@cahoy_transcriptome_2008, we replicate the analysis of @verhaak_integrated_2010

+by using GSVA to transform the gene expression measurements into enrichment

+scores for these four gene sets, without taking the sample subtype grouping

+into account. We start by having a quick glance to the data, which forms part of

+the `r Biocpkg("GSVAdata") package:

+

+```{r}

+library(GSVAdata)

+

+data(gbm_VerhaakEtAl)

+gbm_eset

+head(featureNames(gbm_eset))

+table(gbm_eset$subtype)

+data(brainTxDbSets)

+lengths(brainTxDbSets)

+lapply(brainTxDbSets, head)

+```

+

+GSVA enrichment scores for the gene sets contained in `brainTxDbSets`

+are calculated, in this case using `mx.diff=FALSE`,  as follows:

+

+```{r}

+gbm_es <- gsva(gbm_eset, brainTxDbSets, mx.diff=FALSE, verbose=FALSE)

+```

+

+Figure \@ref(fig:gbmSignature) shows the GSVA enrichment scores obtained for the

+up-regulated gene sets across the samples of the four GBM subtypes. As expected,

+the _neural_ class is associated with the neural gene set and the astrocytic

+gene sets. The _mesenchymal_ subtype is characterized by the expression of

+mesenchymal and microglial markers, thus we expect it to correlate with the

+astroglial gene set. The _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

+_classical_ group correlates highly with the astrocytic gene set. In

+summary, the resulting GSVA enrichment scores recapitulate accurately the

+molecular signatures from @verhaak_integrated_2010.

+

+```{r gbmSignature, height=500, width=700, fig.cap="Heatmap of GSVA scores for cell-type brain signatures from murine models (y-axis) across GBM samples grouped by GBM subtype."}

+library(RColorBrewer)

+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.23,1.21, "Proneural", col="red", cex=1.2)

+text(0.36,1.21, "Neural", col="green", cex=1.2)

+text(0.47,1.21, "Classical", col="blue", cex=1.2)

+text(0.62,1.21, "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)

+```

+

+

+# Parallel calculations

+

+# Frequently asked questions

+

+# Session information {.unnumbered}

+

+Here is the output of `sessionInfo()` on the system on which this document was

+compiled running pandoc `r rmarkdown::pandoc_version()`:

+

+```{r session_info, cache=FALSE}

+sessionInfo()

+```

+

+# References

deleted file mode 100644

@@ -1,960 +0,0 @@

-%\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, echo=FALSE>>=

-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.dif <- gsva(X, gs, mx.diff=TRUE, verbose=FALSE, parallel.sz=1)

-@

-

-\begin{center}

-<<maxvsdif, fig=TRUE, png=TRUE, echo=TRUE, height=5, width=8>>=

-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:

-

-<<results=hide>>=

-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:

-

-<<results=hide>>=

-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:

-

-<<eval=FALSE>>=

-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.

-

-<<figIQR, echo=FALSE, results=hide>>=

-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),

-                           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}.

-

-<<leukemiaVolcano, echo=FALSE, results=hide>>=

-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,