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  . 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)  and (GSA)  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, . Let S denote the indices of the gene set and nS be the size of S. The maxmean statistic S is defined as,
where 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)  , MSigDB  or gene ontology (GO)  . 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  .
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 nS,
where , , , , and . Based on the work in   , an asymptotic distribution of the maxmean statistics under the Lindeberg condition (see  ) 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 f0 by fitting the null gene sets. A special case would be that independently. We can easily compute the parameters:
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 be the density function of the maxmean statistics for all the gene sets. Adopting the two group mixture model   , we specify a similar model that is comprised by a large proportion (p0) of the null density f0 and a small proportion of the non-null density f1,
The null density f0 is assumed to have the form in (2.2), in which the parameters (µ+, µ−, σ+, σ−, ρ) need to be estimated. For identifiability of f0 we further assume the non-null density for some interval A0. Under this assumption, S follows a truncated distribution fT(s) for s A0,
Fitting the maxmean statistics in A0 to (2.4) by maximum likelihood yields . In addition, p0 can be estimated by
In this paper we let A0 be the interval (0, qS) where qS is 90% quantile of the maxmean statistics of all gene sets.
We simulate 2000 gene sets, 5% of them are DE sets and the other 95% are null sets. All sets contain nS = 50 genes. Let C1 and C2 be the sample indices of two conditions and each condition has m = 10 samples. For gene i in sample j, the expression data xij is generated in a hierarchical fashion,
where , for DE sets and δk = 0 for null nets. x0kj is the expression of the hub gene  of set Sk and all genes in the same set are correlated with the hub gene. The DE genes of set Sk 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 f0 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).
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  . 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-
Figure 1. Density estimate by restandardization and empirical method. For null sets, δk = 0, αi = 0 (left) and αi = ±0.2 (right). For DE sets, δk = 1 and αi = ±1. The null proportion is 0.97 (left) and 0.98 (right) estimated by the empirical method.
Figure 2. Estimated null distribution of maxmean statistics of the 4722 GSEA gene sets on the gender dataset. estimated by the empirical method. The point se (red) and sr (blue) are the cutoff values for the empirical and restandardization methods, respectively at FDR = 0.10.
cedure  , 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.
GSA is a representative tool for gene set DE analysis. Compared with the state- of-the-art method GSEA, the maxmean statistic used in GSA is more powerful than the KS statistic in GSEA  . 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 f0 parametrically and all gene sets are compared against f0. The possibility of parametric estimation of f0 is rendered by the large number of gene sets in the hypothesis testing. A similar idea is proposed in  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 correlation 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 A0. In this paper we arbitrarily choose A0 to be the interval (0, qS), where qS is some quantile of the observed maxmean statistics. Determination of qS is a trade-off between variance and bias. With small qS, the bias of the null estimate is small but the variance is large and vice-versa. How to determine an optimal A0 interval requires further exploration.