Estimating the Empirical Null Distribution of Maxmean Statistics in Gene Set Analysis

Gene Set Analysis (GSA) is a framework for testing the association of a set of genes and the outcome, e.g. disease status or treatment group. The method replies on computing a maxmean statistic and estimating the null distribution of the maxmean statistics via a restandardization procedure. In practice, the pre-determined gene sets have stronger intra-correlation than genes across sets. This may result in biases in the estimated null distribution. We derive an asymptotic null distribution of the maxmean statistics based on sparsity assumption. We propose a flexible two group mixture model for the maxmean statistics. The mixture model allows us to estimate the null parameters empirically via maximum likelihood approach. Our empirical method is compared with the restandardization procedure of GSA in simulations. We show that our method is more accurate in null density estimation when the genes are strongly correlated within gene sets.


Introduction
A gene pathway commonly refers to a set of genes that share a particular property, carry out a biological function or lead to a certain product in cells/tissues.Performing differential expression (DE) analysis on such gene sets aggregates the signal of individual genes and potentially increases the power of a hypothesis test.Gene set analysis also provides comprehensive understanding of the biological activities associated with the outcome phenotype and may shed light on treatment of disease.
A variety of tools are available for gene set analysis.These methods can be roughly classified into two broad categories, self-contained and competitive [1].
The self-contained methods test the association between the phenotype and the gene set while ignoring the other genes.Competitive methods overcome this limitation by taking into account other genes when evaluating the association.A common approach is to compute the gene level statistics, and then aggregate them into a gene-set level summary statistic.Of the many competitive methods gene set enrichment analysis (GSEA) [2] and (GSA) [3] are two representative algorithms.In GSEA the Kolmogorov-Smirnov (KS) statistic is employed.GSA improves the power of GSEA by using a more powerful maxmean statistic.First the gene level z statistics are calculated, , 1, , i z i n =  .Let S denote the indices of the gene set and n S be the size of S. The maxmean statistic S is defined as, where {}

I
is the indicator function.
The GSA method estimates the null distribution of S through a restandardization procedure, which is a combination of row (gene) randomization and column (sample label) permutation.It compares the gene set against its permutations and also takes into account the overall distribution of randomly selected null sets.The permutation maintains the correlation structure in the gene set, while row randomization rescales and shifts the permutations to include the competition of genes outside the set.
In practice, the gene sets for analysis are obtained from a public database e.g.
Kyoto Encyclopedia of Genes and Genomes (KEGG) [4], MSigDB [2] or gene ontology (GO) [5].It is reasonable to expect that genes in the same set have stronger correlation than genes across different sets.In such cases there may be a mismatch from the null distribution estimated from restandardization and the true null, which leads to biased inference.
We propose a two group mixture model for the gene set maxmean statistics.
Our model assumes that only a small proportion of the gene sets are truly significantly DE.Based on the mixture model, we apply the maximum likelihood method to estimate the empirical null.The empirical method improves the accuracy of GSA in large scale hypothesis testing.It also reduces the computational burden of the permutation steps in restandardization procedure of GSA.The analysis is demonstrated in simulation studies and a data set from the MSigDB database [2].

Methods
The maxmean statistic in GSA is essentially a maximum statistic of two correlated sample means, S + and S − defined in (1.1), which asymptotically follow bivariate normal distribution for adequately large n S , where ( ) Based on the work in [6] [7], an asymptotic distribution of the maxmean statistics under the Lindeberg condition (see [8]) is where φ and Φ are the probability density function (pdf) and cumulative density function (cdf) of the standard normal distribution.
We can estimate the parameters in f 0 by fitting the null gene sets.A special case would be that We can easily compute the parameters: 0.40 , 0.34 , 0.467.
However, genes in the same set are often correlated.Therefore, the theoretically computed parameters may not match the actual null distribution well.We propose an empirical method to estimate the null distribution of S. Let f be the density function of the maxmean statistics for all the gene sets.Adopting the two group mixture model [9] [10], we specify a similar model that f is com- prised by a large proportion (p 0 ) of the null density f 0 and a small proportion of the non-null density f 1 , ( ) ( ) ( ) 3) The null density f 0 is assumed to have the form in (2.2), in which the parameters (µ + , µ − , σ + , σ − , ρ) need to be estimated.For identifiability of f 0 we further assume the non-null density ( ) for some interval A 0 .Under this assumption, S follows a truncated distribution f T (s) for s ∈ A 0 , Fitting the maxmean statistics in A 0 to (2.4) by maximum likelihood yields 0 f .In addition, p 0 can be estimated by In this paper we let A 0 be the interval (0, q S ) where q S is 90% quantile of the maxmean statistics of all gene sets.

Simulations
We simulate 2000 gene sets, 5% of them are DE sets and the other 95% are null sets.All sets contain n S = 50 genes.Let C 1 and C 2 be the sample indices of two conditions and each condition has m = 10 samples.For gene i in sample j, the expression data x ij is generated in a hierarchical fashion, for DE sets and δ k = 0 for null nets.x 0kj is the expres- sion of the hub gene [11] of set S k and all genes in the same set are correlated with the hub gene.The DE genes of set S k are jointly controlled by δ k and α i .In particular, δ k controls the differential expression of the hub gene between two conditions.The parameter α i controls the inter-correlation within the set.When α i =0 genes are independent within gene sets.For null sets we let α i = 0 (independent) or ±0.2 (correlated).For DE sets α i = ±1 so that the correlation is stronger than the null sets.Figure 1 shows the f 0 density estimate by our empirical method and the restandardization procedure in GSA.
When genes within the gene sets are independent, the null distribution obtained by the two methods are very close (Figure 1 left).As 95% of the gene sets are null sets, the estimated distributions match the majority of the sets.When genes are correlated, the null distribution obtained by GSA shows a mismatch from the null gene sets, while our empirical method maintains a good fit to the data (Figure 1 right).

Application
We evaluate our empirical method and GSA with the 4722 curated gene sets in MSigDB database using the gender dataset included in the GSEA package [2].
The minimum size of the gene sets is 5, the maximum size is 1972 and median size is 39.The dataset contains transcriptional profiles of 15,056 genes in 17 female and 15 male lymphoblastoid cell line samples.There are 302 gene sets with unadjusted p-value < 0.05 by our empirical method and 507 gene sets by GSA.
Controlling the false discovery rate (FDR) at 0.1 with Benjamini-Hochberg pro-  cedure [12], the empirical method identifies 39 significant sets while the restandardization identifies 8 significant sets.Figure 2 shows the null distribution obtained by the empirical method and GSA.

Discussion
GSA is a representative tool for gene set DE analysis.Compared with the stateof-the-art method GSEA, the maxmean statistic used in GSA is more powerful than the KS statistic in GSEA [3].GSA also establishes a restandardization algorithm to assess the maxmean statistic.Due to the correlation of the genes in pre-defined gene sets or pathways, the null distribution obtained by the restandardization procedure in GSA may not match the true null distribution well.
In studying the GSA method, we propose a new method to estimate the null distribution of the gene set maxmean statistics.Unlike the permutation test of GSA in which every gene set is compared against its own permutations, the fundamental difference of our method is that it estimates an overall null distribution f 0 parametrically and all gene sets are compared against f 0 .The possibility of parametric estimation of f 0 is rendered by the large number of gene sets in the hypothesis testing.A similar idea is proposed in [9] for large-scale test of DE genes.
We extend the idea to the test of gene sets when a large number of sets are available.
Our method is based on the sparsity assumption that only a small proportion of the gene sets are truly DE.Further, we adopt a two group parametric mixture distribution to model the gene set maxmean statistics.The parameters of the null distribution is estimated under the sparsity assumption.We show that our new method provides more accurate estimation of the null distribution.It also avoids the computational intensity in the permutation steps of GSA.
In the simulations, we compare the two methods under independence and correlation.When the intra-set correlation is greater than the cross-set correlation, the GSA shows a mismatch to the true null.The reason is that the randomization step samples genes in the entire dataset, which has a different correla-tion structure from gene sets.As a result the mean and standard deviation obtained in the randomization step is inaccurate.
In application to the gender data set, our method has fewer gene sets with unadjusted p-value < 0.05, but identifies more gene sets after controlling for FDR.A reason is that in GSA, gene-set p-values are limited by the number of permutations.Due to the large number of genes and gene sets, with 10,000 permutations, GSA takes 64 minutes on an Intel i7 3770K 3.5 GHz CPU.In contrast, our empirical method takes less than 1 minute.
An important aspect of our empirical estimation method is specifying an appropriate range for the zero assumption region A 0 .In this paper we arbitrarily choose A 0 to be the interval (0, q S ), where q S is some quantile of the observed maxmean statistics.Determination of q S is a trade-off between variance and bias.With small q S , the bias of the null estimate is small but the variance is large and vice-versa.How to determine an optimal A 0 interval requires further exploration.

Figure 2 .
Figure 2.Estimated null distribution of maxmean statistics of the 4722 GSEA gene sets on the gender dataset.0 ˆ0.96 p = estimated by the empirical method.The point s e (red) and s r (blue) are the cutoff values for the empirical and restandardization methods, respectively at FDR = 0.10.