GSAR: Bioconductor package for Gene Set analysis in R

Background Gene set analysis (in a form of functionally related genes or pathways) has become the method of choice for analyzing omics data in general and gene expression data in particular. There are many statistical methods that either summarize gene-level statistics for a gene set or apply a multivariate statistic that accounts for intergene correlations. Most available methods detect complex departures from the null hypothesis but lack the ability to identify the specific alternative hypothesis that rejects the null. Results GSAR (Gene Set Analysis in R) is an open-source R/Bioconductor software package for gene set analysis (GSA). It implements self-contained multivariate non-parametric statistical methods testing a complex null hypothesis against specific alternatives, such as differences in mean (shift), variance (scale), or net correlation structure. The package also provides a graphical visualization tool, based on the union of two minimum spanning trees, for correlation networks to examine the change in the correlation structures of a gene set between two conditions and highlight influential genes (hubs). Conclusions Package GSAR provides a set of multivariate non-parametric statistical methods that test a complex null hypothesis against specific alternatives. The methods in package GSAR are applicable to any type of omics data that can be represented in a matrix format. The package, with detailed instructions and examples, is freely available under the GPL (> = 2) license from the Bioconductor web site. Electronic supplementary material The online version of this article (doi:10.1186/s12859-017-1482-6) contains supplementary material, which is available to authorized users.


Introduction
This vignette provides an overview of the R package GSAR which provides a set of multivariate statistical tests for self-contained gene set analysis (GSA). GSAR consists of two-sample multivariate nonparametric statistical methods testing a null hypothesis against specific alternative hypotheses, such as differences in mean (shift), variance (scale) or correlation structure. It also offers a graphical visualization tool for the correlation networks obtained from expression data to examine the change in the net correlation structure of a gene set between two conditions based on the minimum spanning trees. The same tool also works for protein-protein Interaction (PPI) networks to highlight the most influential proteins. The package implements the methods proposed in [1, 2,3] which were thoroughly tested using simulated and microarray datasets in [1] and [2]. Figure 1 shows the outline of the package. While the test functions in the package analyze one gene set at a time, the wrapper function TestGeneSets receives a list of gene sets and invokes a specific test function for all the gene sets in a sequential manner.
The methods in the package can also be applied to RNA-seq count data given that proper normalization which accounts for both the within-sample differences (gene lengths) and between-samples differences (library sizes) is used. However, for RNA-seq data the variance is a function of mean expression and applying the variance tests (RKStest, RMDtest, AggrFtest) to RNA-seq data highly depends on the normalization method. This specific topic is not well-studied and characterized yet. This vignette begins with a brief overview of the theoretical concepts behind the statistical methods, and then provides a number of fully worked case studies, from microarrays to RNA-seq count data. MD statistic (discrete Normal-like distribution) Fisher-combined F-statistic (continuous Chi-squarelike distribution) W GSNCA statistic (continuous Normallike distribution) Figure 1: GSAR package outline. The inputs for the statistical tests can be (1) the matrix of gene expression for a single gene set in the form of normalized microarray or RNA-seq data and a vector of labels indicating to which condition each sample belongs; or (2) the matrix of gene expression for all genes, a vector of labels indicating to which condition each sample belongs, and a list of gene sets. Each test returns P-value and, optionally, the test statistic of observed data, test statistic for all permutations, and other optional outputs. Some functions produce graph plots.
Many methodologies for testing the differential expression of gene sets have been suggested and are collectively named gene set analysis (GSA). GSA can be either competitive or self-contained. Competitive approaches compare a gene set against its complement which contains all genes excluding the genes in the set, and self-contained approaches compare whether a gene set is differentially expressed (DE) between two phenotypes (condition). Some competitive approaches are influenced by the genomic coverage and the filtering of the data and can increase their power by the addition of unrelated data and even noise [4] while others can be influenced by the proportion of up-regulated and down-regulated genes in a gene set between two conditions [5]. Due to these problems, package GSAR focuses on self-contained methods only. The possibility to formulate different statistical hypotheses by using different test statistics with self-contained approaches enables the formulation and exploration of different biological hypotheses [2]. For GSA, testing hypotheses other than the equality of the mean expression vectors remains underexplored.
Package GSAR provides a set of methods to test a null hypothesis against specific alternatives, such as differential distribution (function WWtest), shift or mean (functions KStest and MDtest), scale or variance (functions RKStest, RMDtest, and AggrFtest) or net correlation structure (function GSNCAtest).
Most of the statistical methods available in package GSAR (all except GSNCAtest and AggrFtest) are network or graph-based. GSAR handles graphs using the rich functionality of the igraph class from package igraph [6]. GSAR invokes some functions from package igraph in its methods implementation and uses the plot method for class igraph for visualizing the generated graphs. In what follows we introduce the following notations. Consider two different biological phenotypes, with n 1 samples of measurements for the first and n 2 samples of the same measurements for the second. Each sample is a p-dimensional vector of measurements of p genes (constituting a single gene set). Hence, the data for the first phenotype consists of a p × n 1 matrix and for the second phenotype consists of a p × n 2 matrix, where rows are genes and columns are samples. The samples of the first and second phenotypes are respectively represented by two random vectors X and Y . Let X and Y be independent and identically distributed with the distribution functions F x , F y , p-dimensional mean vectorsX andȲ , and p × p covariance matrices S x and S y .
2 Minimum spanning trees

First MST
The pooled multivariate (p-dimensional) observations X and Y can be represented by an edge-weighted graph G(V, E) where V is the set of vertices in the graph. Each vertex in the sample network corresponds to one observation (sample) and E is the set of edges connecting pairs of vertices. The complete graph of X and Y has N = n 1 + n 2 vertices and N (N − 1)/2 edges. The weights of the edges are estimated by the Euclidean distances between pairs of observations (samples) in R p .
The minimum spanning tree (MST) is defined as the acyclic subset T 1 ⊆ E that connects all vertices in V and whose total length i,j∈T1 d(v i , v j ) is minimal. Each vertex in the graph corresponds to a p-dimensional observation from X or Y . The MST provides a way of ranking the multivariate observations by giving them ranks according to the positions of their corresponding vertices in the MST. The purpose of this ranking is to obtain the strong relationship between observations differences in ranks and their distances in R p . The ranking algorithm can be designed specifically to confine a particular alternative hypothesis more detection power [2]. Five tests in package GSAR are based on MST: WWtest, KStest, MDtest, RKStest, and RMDtest.
The following example generates a feature set of 20 features and 40 observations using the random multivariate normal data generator from package MASS, creates a graph object from the data and obtain its MST using functions from package igraph.

MST2 for correlation and PPI networks
The second MST is defined as the MST of the reduced graph G(V, E − T 1 ). We denote the union of the first and second MSTs by MST2. Each vertex in the MST2 has a minimum degree of 2 if all the edges between vertices are considered (full network).
The correlation (coexpression) network is defined as the edge-weighted graph G(V, E) where V is the set of vertices in the graph with each vertex corresponding to one feature (gene) in the gene set and E is the set of edges connecting pairs of vertices with weights estimated by some correlation distance measure. The correlation distance here is defined by d ij = 1 − |r ij | where d ij and r ij are respectively the correlation distance and correlation coefficient between genes i and j [1]. The MST2 of the correlation network gives the minimal set of essential links (interactions) among genes, which we interpret as a network of functional interactions. A gene that is highly correlated with most of the other genes in the gene set tends to occupy a central position and has a relatively high degree in the MST2 because the shortest paths connecting the vertices of the first and second MSTs tend to pass through this gene. In contrast, a gene with low intergene correlations most likely occupies a non-central position in the MST2 and has a degree of 2. This property of the MST2 makes it a valuable graphical visualization tool to examine the full correlation network obtained from gene expression data by highlighting the most influential genes. As an example, the MST2 of the dataset generated in the previous example can be found as follows > ## The input of findMST2 must be a matrix with rows and columns > ## respectively corresponding to genes and columns. > ## Therefore, dataset must be transposed first. > dataset <-aperm(dataset, c(2,1)) > MST2 <-findMST2(dataset) The protein-protein Interaction (PPI) network is defined as the binary graph G(V, E) where V is the set of vertices in the graph with each vertex corresponding to one protein and E is the set of edges connecting pairs of vertices where e ij = 1 if interaction exists between proteins i and j and e ij = 0 otherwise. The MST2 of the PPI network gives the minimal set of essential interactions among proteins. It was shown in [7] that function plotMST2.PPI reveals fine network structure with highly connected proteins occupying central positions in clusters.

Statistical methods
Package GSAR provides seven statistical methods that test five specific alternative hypotheses against the null (see Figure 1), two functions to find the MST2 of correlation and PPI networks, one wrapper function to plot the MST2 of correlation networks obtained from gene expression data under two conditions, and one wrapper function to facilitate performing a specific statistical method for a list of gene sets.

Wald-Wolfowitz test
The Wald-Wolfowitz (WW) tests the null hypothesis H 0 : F x = F y against the alternative H 1 : F x = F y . When p = 1, the univariate WW test begins by sorting the observations from two phenotypes in ascending order and labeling each observation by its phenotype. Then, the number of runs (R) is calculated where R is a consecutive sequence of identical labels. The test statistic is a function of the number of runs and is asymptotically normally distributed.
The multivariate generalization (p > 1) suggested in [3] is based on the MST. Similar to the univariate case, in the multivariate generalization of WW test, all edges in the MST connecting two vertices (observations) with different labels are removed and the number of the remaining disjoint trees (R) is calculated [3]. The test statistic is the standardized number of subtrees The null distribution of W is obtained by permuting the observation labels for a large number of times (nperm) and calculating W for each time. The null distribution is asymptotically normal. P -value is calculated as where W k is the test statistic for permutation k, W obs is the observed test statistic, and I is the indicator function. Function WWtest performs this test.

Kolmogorov-Smirnov tests
When p = 1, the univariate Kolmogorov-Smirnov (KS) test begins by sorting the observations from two phenotypes in ascending order. Then observations are ranked and the quantity The test statistic is the maximal absolute difference between the Cumulative Distribution Functions (CDFs) of the ranks of samples from X and Y, D = n1n2 n1+n2 max|d i |. The null distribution of D is obtained by permuting the observation labels for a large number of times (nperm) and calculating D for each time. The KS statistic asymptotically follows the Smirnov distribution and tests a one-sided hypothesis. P -value is calculated as where D k is the test statistic for permutation k, D obs is the observed test statistic, and I is the indicator function.
The ranking scheme can be designed to confine a specific alternative hypothesis more power. Two possibilities are available: First, if the null hypothesis H 0 :μ X =μ Y is tested against the alternative H 1 :μ X =μ Y , the MST is rooted at a node with the largest geodesic distance and the rest of the nodes are ranked according to the high directed preorder (HDP) traversal of the tree [3], which can be found using function HDP.ranking. Function KStest performs this specific test. Second, if the null hypothesis H 0 :σ X =σ Y is tested against the alternative hypothesis H 1 :σ X =σ Y , the MST is rooted at the node of smallest geodesic distance (centroid) and nodes with largest depths from the root are assigned higher ranks. Hence, ranks are increasing radially from the root of the MST. The radial ranking of vertices in a tree can be found using function radial.ranking. Function RKStest performs this specific test.
The MST found in the previous example is shown in Figure 2 were vertices from group 1 are in green color and vertices from group 2 are in yellow color. Ranking vertices in the graph according to the HDP and radial traversal of the MST can be done respectively using functions HDP.ranking and radial.ranking

Mean deviation tests
The mean deviation (MD) statistic calculates the average deviation between the Cumulative Distribution Functions (CDFs) of the ranks of samples from X and Y. The MD statistic for a gene set of size p is defined as r j is the rank of sample j in the MST and the exponent α is set to 0.25 to give the ranks a modest weight. The null distribution of D is obtained by permuting sample labels for a large number of times (nperm) and calculating D for each time. The MD statistic has asymptotically a normal distribution and tests a two-sided hypothesis. P -value is calculated as where D k is the test statistic for permutation k, D obs is the observed test statistic, and I is the indicator function. Similar to KS statistic, combining this statistic with appropriate ranking schemes, allows the test of specific alternative hypotheses, namely diffferential mean (mean deviation test function MDtest) and differential variance (radial mean deviation function RMDtest).

Aggregated F-test
The univariate F-test is used to detect differential variance in individual genes forming a gene set of size p, and the individual P-values, p i , 1 ≤ i ≤ p, are then aggregated using Fisher probability combining method to obtain a score statistic (T ) for the gene set Significance is estimated by permuting sample labels and calculating T for a large number of times (nperm). P-value is calculated as where T k is the test statistic for permutation k, T obs is the observed test statistic, and I is the indicator function.

Gene sets net correlations analysis 3.5.1 Method
Gene sets net correlations analysis (GSNCA) is a two-sample nonparametric multivariate differential coexpression test that accounts for the correlation structure between features (genes). For a gene set of size p, the test assigns weight factors w i , 1 ≤ i ≤ p, to genes under one condition and adjust these weights simultaneously such that equality is achieved between each genes's weight and the sum of its weighted absolute correlations (r ij ) with other genes in a gene set of p genes The problem is solved as an eigenvector problem with a unique solution which is the eigenvector corresponding to the largest eigenvalue of the genes' correlation matrix (see [1] for details).
The test statistic w GSN CA is given by the first norm between the scaled weight vectors w (1) and w (2) (each vector is multiplied by its norm) between two conditions This test statistic tests the null hypothesis H 0 : w GSN CA = 0 against the alternative H 1 : w GSN CA = 0. The performance of this test was thoroughly examind in [1]. P -value is calculated in exactly the same way as before for the WW and KS tests. The values in the scaled weight vectors w (1) and w (2) roughly fall in the range [0.5, 1.5], with high values indicating genes that are highly correlated with other genes in the same gene set.

The problem of zero standard deviation
In special cases some features in a set may have constant or nearly constant levels across the samples in one or both conditions. Such situation almost never encountered in microarray data, but may arises for RNA-seq count data where a gene may have zero counts under many samples if the gene is not expressed. This results in a zero or a tiny standard deviation. Such case produces an error in command cor used to compute the correlations between features.
To avoid this situation, standard deviations are checked in advance (default behaviour) and if any is found below a specified minimum limit (default is 1e-3), the execution stops and an error message is returned indicating the the number of feature causing the problem (if only one the index of that feature is given too). To perform the GSNCA for count data, the features causing the problem must be excluded from the set.
If a feature has nearly a constant level for some (but not all) samples under both conditions, permuting sample labels may group together such samples under one condition by chance and hence produce a standard deviation smaller than the minimum limit. To allow the test to skip such permutations without causing excessive delay, an upper limit for the number of allowed skips can be set (default is 10). If the upper limit is exceeded, an error message is returned.
If the user is certain that the tested feature sets contain no feature with nearly zero standard deviation (such as the case with filtered microarray data), the checking step for tiny standard deviations can be skipped in order to reduce the execution time.
4 Notes on handling RNA-seq count data RNA-seq data consists of integer counts usually represented by the discrete Poisson or Negative Binomial distributions. Therefore, tests designed for microarray data (which follows the continuous normal distribution) can not be applied directly to RNA-seq data. The nonparametric tests presented in package GSAR need no prior distributional assumptions and can be applied to RNA-seq counts given that proper normalization is used. The normalization should accounts for the between-samples differences (library size or sequencing depth) and within-sample differences (mainly gene length). The reads per kilobase per million (RPKM) is such normalization. However, due to the lack of thorough performance studies, two points must be declared: The variance of both the Poisson and negative Bionomial distributions, used to model RNA-seq count data is a function of their mean. Therefore, applying the variance tests to RNA-seq data highly depends on the normalization method. This specific topic is not well-studied and characterized yet. RNA-seq datasets often have many zero counts, therefore, the problem of having at least one gene with zero standard deviations in a gene set is frequent and prevent calculating the correlation coefficients necessary to perform the GSNCA. One possible solution is to discard any genes that may have zero or tiny standard deviation and apply GSNCA to the remaining genes in the gene set.

Case studies
This Section illustrates the typical procedure for applying the methods available in package GSAR to perform GSA. Two microarray and one RNA-seq datasets are used.  [8,9]. Transcriptional profiles obtained from microarrays of platform hgu95av2 are available from the Broad Institute's website (http://www.broadinstitute.org/gsea/datasets.jsp).

Filtering and normalization
A preprocessed version of p53 dataset is available in package GSAR as a matrix object. The p53 dataset was dowloaded from the Broad Institute's website. Probe level intensities were quantile normalized and transformed to the log scale using log 2 (1 + intensity). Probes originally had Affymetrix identifiers which are mapped to unique gene symbol identifiers. Probes without mapping to entrez and gene symbol identifiers were discarded. Probes with duplicate intensities were assessed and the probe with the largest absolute value of t-statistic between WT and MUT conditions was selected as the gene match. Finally, genes were assigned gene symbol identifiers and columns were assigned names indicating weither they belong to WT or MUT group. The columns were sorted such that the first 17 columns are WT samples and the next 33 columns are the MUT samples. This processed version of the p53 dataset was used in the analysis presented in [1].

GSA
GSA is performed on selected C2 curated gene sets (pathways) of the molecular signatures database (MSigDB) 3.0 [10]. This list of gene sets is available in package GSVAdata. We start by loading the required data > library(GSVAdata) > data(p53DataSet) > data(c2BroadSets) c2BroadSets is an object of class GeneSetCollection supported by package GSEABase. The genes in the c2BroadSets object have entrez identifiers. Package org.Hs.eg.db is used to convert the entrez identifiers to gene symbol identifiers. Genes without unique mapping to gene symbol identifiers or that do not exist in the p53 dataset are discarded from the C2 pathways. This insures proper indexing of genes in the dataset by the gene names in each C2 pathway. Finally, we keep only pathways with 10 ≤ p ≤ 500 where p is the number of genes remaining in the pathways after filtering steps.
> library(org.Hs.eg.db) > library(GSEABase) > C2 <-as.list(geneIds(c2BroadSets)) > len <-length(C2) > genes.entrez <-unique(unlist(C2)) > genes.symbol <-array("",c(length(genes.entrez),1)) > x <-org.Hs.egSYMBOL > mapped_genes <-mappedkeys( The questions addressed by these tests were the identification of gene sets expressed with different distributions, means, variances or correlation structure between two conditions. At a significance level 0.05, the targeted pathway shows a statistical evidence of being differentially coexpressed only. The MST2s of the correlation network for WT and MUT groups are shown in Figure 3, generated by function plotMST2.pathway  (1,17), rep(2,33)), name="LU_TUMOR_VASCULATURE_UP", + legend.size=1.2, leg.x=-1.2, leg.y=2, + label.size=1, label.dist=0.8, cor.method="pearson") The targeted pathway comprises genes over-expressed in ovarian cancer endothelium [11]. Gene TNFAIP6 (tumor necrosis factor, α-induced protein 6) identified by GSNCA as a hub gene for WT group and visualized using MST2 (Figure 3, left panel) was found 29.1-fold over-expressed in tumor endothelium in the original study and was suggested to be specific for ovarian cancer vasculature. This indicates that gene TNFAIP6 can be an important regulator of ovarian cancer, and identifying it as a hub by GSNCA enhances the original observation. When p53 is mutated (Figure 3, right panel) the hub gene is VCAN, containing p53 binding site and its expression is highly correlated with p53 dosage [12]. Therefore, both hub genes provide adequate information about the underlying biological processes. If testing multiple or all the gene sets in the c2.pathways list is desired, the wrapper function TestGeneSets can be used. For example, the following code tests the first 3 gene sets in list c2.pathways, that have a minimum of 10 and a maximum of 100 genes, using the GSNCA method This dataset consists of microarrays (platform hgu95av2) from 128 different individuals with acute lymphoblastic leukemia (ALL). There are 95 samples with B-cell ALL [13] and 33 with T-cell ALL [14]. We consider B-cell type only and compare tumors carrying the BCR/ABL mutations (37 samples) to those with no cytogenetic abnormalities (42 samples). The Bioconductor package ALL provides the ALL dataset with samples normalized using the robust multiarray analysis (RMA) procedure [15].

Filtering and normalization
Affymetrix probe identifiers without mapping to entrez and gene symbol identifiers were discarded. Affymetrix probe identifiers were mapped to unique gene symbol identifiers and intensities of probes mapping to the same gene symbol identifier were assessed and the probe with the largest absolute value of t-statistic between normal (NEG) and mutation (MUT) conditions was selected as the gene match.

Selected gene set
Lets examine the C2 pathway KEGG CHRONIC MYELOID LEUKEMIA, knowm to be specifically associated with the BCR/ABL mutation. This pathway has many BCR/ABL-related genes and hence expected to show difference between NEG and MUT conditions. To ensure proper indexing, the lists of genes in C2 pathways should consists only of genes available in the filtered ALL dataset. Therefore, the same steps taken to filter the C2 pathways with the p53 dataset are repeated for the ALL dataset.

Introduction
The Pickrell dataset of sequenced cDNA libraries generated from 69 lymphoblastoid cell lines derived from unrelated Yoruban Nigerian individuals (YRI) is part of the HapMap project. The original experimental data was published by [16]. Package tweeDEseqCountData provides the table of counts for this dataset in the expression set object pickrell.eset. This table of counts corresponds to the one in the ReCount repository available at http:// bowtie-bio.sourceforge.net/recount. Details on the pre-processing steps to obtain this table of counts from the raw reads are provided by [17].
Package tweeDEseqCountData provides annotation data for the human genes forming the table in pickrell.eset as a data frame object annotEnsembl63. tweeDEseqCountData also provides two lists of genes (gene sets) with documented sex-specific expression and occurring within the set of genes that form the table of counts in pickrell.eset. The first is a set of genes that are located on the male-specific region of chromosome Y, and therefore are over-expressed in males (msYgenes). The second is a set of genes, that are escaping X-chromosome inactivation, and therefore are overexpressed in females (XiEgenes). These two sets are useful in serving as true positives when GSA is conducted between males and females to detect gene sets that are differentially expressed.

Filtering and normalization
Any transcript without entrez identifier or gene length information is discarded. To consider only expressed genes in the analysis, genes with an average count per million (cpm) less than 0.1 are discarded. The gene length information is used to perform the RPKM normalization. Finally, the RPKM-normalized expression is transformed to the logarithm scale. RPKM as well as a few other normalizations were used with the Pickrell datasets in [18] to perform GSA and the study found no significant differences between different normalizations for the same test statistic.

Testing selected pathways
Any gene in msYgenes, XiEgenes, or Xigenes but not found in the filtered dataset is discarded. Then, the remaining gender-related genes in msYgenes and XiEgenes are combined into one gene set (XYgenes [1] 147 Notice that two genes (ENSG00000183878 and ENSG00000198692) in XYpathway had zero standard deviations for female samples and were filtered out from XYpathway. These two genes are Y-liked genes and expected to have many zero counts for female samples. Although this filtering step increases the chances of success in performing GSNCA, the existence of many zero counts dispersed over many samples for one or more genes may still cause a problem when the sample permutation process groups many zero counts under one condition. The parameter max.skip in function GSNCAtest allows some tolerance by assigning the maximum number of skipped permutations allowed to avoid the few ones causing the problem. This solution may work or fail depending on the proportion of zero counts in the data. For example, assigning max.skip to 100 or more solved the problem for XYpathway, but it did not for Xipathway. We advise to perform gene filtering based on zero counts prior to trying the GSNCA for count data.