s2fieldsofstudy: Biology Chemistry

text:
A novel multikinase inhibitor R8 exhibits potent inhibition on cancer cells through both apoptosis and autophagic cell death

Chemotherapy is an important treatment for cancer patients, especially for those with unresectable lesions. Targeted therapy of cancer by specific inhibition of aberrant tyrosine kinase activities in cancer cells with chemically synthesized tyrosine kinase inhibitors (TKIs), shows better responses while less side effects than traditional chemotherapeutic drugs. It is common that cancer cells often exhibit deregulation of several tyrosine kinases simultaneously, multikinase TKIs (MKIs) therefore have greater advantages over single-target TKIs. Currently more MKIs are under developing for better efficacy for different types of cancer. In the present work, we evaluated the in vitro therapeutic potential of a novel MKI, namely R8, with comparison to the clinically available MKI Sunitinib. Results showed that R8 has stronger inhibition on six different types of cancer cell lines with lower IC50 than Sunitinib does. Cell cycle analysis showed that R8 induced significant G0/G1 arrest phase of lung cancer A549 and NCI-H226 cells. The inhibition was also confirmed by colony formation and migration assays in both lung cancer cell lines in a dose-dependent manner. R8 could significantly inhibit the phosphorylation of multiple receptor tyrosine kinases (RTKs) included PDGFRβ, VEGFR2, EGFR and C-Kit, leading to the down-regulation of PI3K-Akt-mTOR signaling. Further analysis revealed that R8 treatment induced more significant apoptosis than Sunitinib did, which might be the consequence of the autophagic cell death. In conclusion, this work suggested R8 to be a promising novel anticancer MKI, and provided the basis for further in vivo investigation on its potential in treatment of lung cancer.

INTRODUCTION
Cancer is one of the leading causes of death worldwide, with an estimated 14.1 million new cases and 8.2 million deaths occurred in 2012 worldwide [1]. Chemotherapy is an important way of treatment for cancer patients, especially for those with unresectable lesions. Targeted therapy of cancer by specific inhibition of aberrantly activated receptor tyrosine kinase (RTK) activities in cancer cells with chemically synthesized small molecule compounds, tyrosine kinase inhibitors (TKIs), shows better responses while less side effects than traditional chemotherapeutic drugs. It is common that cancer cells often exhibit simultaneous deregulation of several tyrosine kinases, such as epithelial growth factor receptor (EGFR), platelet-derived growth factor receptor (PDGFR), vascular endothelial growth factor receptor (VEGFR), etc., which lead to the malignant characteristics www.impactjournals.com/oncotarget/ Oncotarget, 2017, Vol. 8, (No. 50), pp: 87209-87220 Research Paper of tumor such as uncontrolled proliferation, invasiveness and angiogenesis, multi-target TKIs (or multikinase inhibitors, MKIs) therefore have great advantages over single-target TKIs [2].
Although many TKIs/MKIs have been widely used clinically in the treatment of various types of cancer, they do not always give satisfactory outcomes, even for individuals with definite hyper-activated target kinases. Primary resistance and acquired resistance are commonly observed because of the fact that, each type of cancer or even each individual cancer cell may show different aberrantly activated tyrosine kinases due to heterogenicity during the development and progression of the disease, and the TKIs/MKIs may also have different binding and inhibition selectivities on different kinases (and even same kinase with different mutations) [3][4][5]. Therefore more multi-target TKIs are developed or under developing for better efficacies on different individuals and to overcome the primary or acquired resistance [6][7][8].
In the present work, we evaluated the in vitro inhibitory effect of a novel multikinase inhibitor, namely R8 (Figure 1), on several different types cancer cells. We further focused on the mechanism for its inhibition on lung cancer cells specifically, which showed distinct characteristics compared to Sunitinib, a clinically available MKI which has been investigated extensively in many types of cancer including lung cancer, either used alone or in combination with other therapeutics [9][10][11].

R8 significantly inhibited proliferation of nonsmall lung cancer cell line in a dose-dependent manner
In order to evaluate efficacy of R8 on inhibiting cell proliferation, we selected clinically available MKI Sunitinib as a control, which was reported to be successful in inhibiting proliferation of variety kinds of tumors. We selected six cell lines, list as follows: human gastric carcinoma cells SGC7901, non-small lung cancer cell line A549 and NCI-H226, laryngeal carcinoma cell line HEp-2, renal adenocarcinoma cell line 786-0, colorectal cancer cell line SW620. As indicated by the cell viability curves (Figure 2A- Figure 2F)   Higher G0/G1 phase cell cycle arrest was induced after A549 and NCI-H226 cell lines treated with R8 for 24 h than Sunitinb. Data were presented as means ± SD of three independent experiments ( * P<0.05, ** P<0.01, *** P<0.001). www.impactjournals.com/oncotarget Furthermore, the long-term effect of R8 on cell survival was determined by colony formation assay. The results further confirmed that R8 was more significantly in inhibiting the proliferation of both A549 ( Figure 2G- Figure 2I) and NCI-H226 ( Figure 2J- Figure 2L) than Sunitinb.

R8 induced G0/G1-phase cell cycle arrest and inhibited cell migration in non-small lung cancer cells
To identify the mechanism underlying R8-mediated inhibition of proliferation, cell cycle analysis was performed. R8 caused enhanced accumulation of cells in G1 phase and decrease in S phase than Sunitinib in A549 ( Figure 2M) and NCI-H226 ( Figure 2N). Significantly induction of G0/G1 phase arrest may be one of reasons to answer why R8 was more competent in inhibiting cell proliferation. Since migration and angiogenesis are essential for tumor metastasis, so next we investigated whether R8 could inhibit cell migration ability. As shown in Figure 3, R8 significantly inhibited NCI-H226 ( Figure  3A, Figure 3B) and A549 ( Figure 3C) migration in a dose-dependent manner than Sunitinib. But when we evaluated tube formation of HUVECs treated with Sunitinib or R8, R8 was less efficient in inhibiting VEGF pathway than Sunitinib ( Figure 3D, Figure 3E).

R8 inhibited growth factor stimulated receptor phosphorylation
Above data indicated R8 could inhibit A549 and NCI-H226 cells proliferation and migration, but the IC 50 values in two cell lines were different, NCI-H226 cell line seemed more resistant to drug treatment. So next we compared the expression levels of PDGFRβ, VEGFR2, EGFR, C-kit between A549 and NCI-H226. Undoubtedly, NCI-H226 expressed high level of totally four receptors, especially for EGFR and PDGFRβ, whereas A549 only with medium level of EGFR and low level of C-Kit expression ( Figure 4A). So next, we evaluated the phosphorylated status of tyrosine kinase receptors in A549 and NCI-H226 after exposed to Sunitinib or R8 for 1 h, respectively. Sunitinib treatment showed slight decrease in p-EGFR but no change with p-C-Kit, whereas R8 (10 μM) treated A549 caused an up-regulation in p-EGFR. P-AKT and p-ERK were activated after stimulation of 10% FBS, and Sunitinib inhibited phosphorylation of p-ERK in a dose-dependent manner whereas no change with p-AKT. Although R8 caused up-regulation of p-EGFR, dosedependent inhibition of p-AKT, total AKT, p-ERK, total ERK were observed in R8 treated A549 cell line ( Figure  4B).
For NCI-H226, which expressed abundant receptors, 10% FBS stimulation for 10 min caused significant upregulation of four receptors. Sunitinib and R8 inhibited p-PDGFRβ, p-VEGFR2, p-EGFR, p-C-Kit expression in a dose-dependent manner. Although Sunitinib and R8 downregulated p-EGFR and p-C-Kit, but was not so effective than p-PDGFRβ, p-VEGFR2. P-Akt and total AKT were decreased in R8 (20 μM) treated NCI-H226 cells, but only a slight decrease of p-AKT in Sunitinib (20 μM) treated groups ( Figure 4C).

R8 induced autophagy and apoptosis
Since treated with Sunitinib and R8 caused a decline of phosphorylated receptors which may induce cell starvation, accompanied with down-regulated p-AKT, both of which may result in autophagy.
So next we analyzed autophagy and apoptosis related proteins by Western blotting. LC3-II was a good indicator of autophagosome formation. For A549, R8 (5 μM, 10 μM) induced significant increase of LC3-II than Sunitinib in a dose-dependent manner. Apoptosis related proteins, such as cleaved-caspase-3 and cleaved PARP, were observed after treated with Sunitinib or R8 for 48 h, but no significant up-regulation of Bax, a pro-apoptotic molecule ( Figure 5A). For NCI-H226, the expression of LC3-II was up-regulated in a dose-dependent manner both in sunitinb and R8 treated group, whereas R8 showed more potent in accumulating of LC3-II. Significant cleavage of caspase-3 and PARP, up-regulation of Bax were observed in R8 treated cells than Sunitinib ( Figure 5B).
Since there were many vacuole-like bubbles in cytoplasm after NCI-H226 treated with R8 for 48 h compared with controls, so we performed LC3 immunofluorescence to distinguish cytosolic LC3-I and autophagosome-associated LC3-II. For NCI-H226 exposed to 20 μM R8, red-punctate LC3-II was observed in the cytoplasm, which indicated the formation of autophagosome ( Figure 5C).

R8 induced type II programmed cell death
In both A549 and NCI-H226 cell lines, R8 induced significant autophagy and apoptosis than Sunitinib. But the role of autophagy induced by Sunitinib or R8 was not clear. So we first incubated cells with autophagy inhibitor 3-Methyladenine (3-MA) or chloroquine phosphate (CQ) before drug treatment.
For A549, inhibition of autophagy by 3-MA or CQ resulted in increasing of cleaved caspase-3 and PARP, which indicated autophagy induced by Sunitinib in A549 was a self-protected autophagy. Once inhibited autophagy, apoptosis would increase. On the other hand, for R8 treated cells, 3-MA pre-treatment resulted in decreasing of Beclin-1, LC3-II, cleaved caspase-3 and PARP, which indicated that the autophagy induced by R8 was deathleading ( Figure 6A). Once inhibited autophagy, cell would survive. CQ pre-treated A549 only showed slight function in inhibiting cell-death autophagy, which indicated class III-PI3K was more important in R8-induced autophagy.
So in NCI-H226, CQ was more potent in inhibiting R8 induced autophagy, which indicated fusion with lysosome was more important in R8-induced autophagy of NCI-H226 cell line.
Besides, cell viability was up-regulated after A549 ( Figure 6C) and NCI-H226 ( Figure 6D) treated with combination of R8 and 3-MA, but not for the Sunitinib treated cells.
In conclusion, Sunitinib induced a self-protected autophagy, protected cells from apoptosis; whereas R8 induced autophagic cell death, which contributed to apoptosis.

DISCUSSION
In this study, we observed the inhibition of a novel MKI R8 on cancer cells in vitro, with comparison to Sunitinib, one of the widely used MKIs in clinic, and explored the possible mechanism by which R8 took effect. As a class III/V multikinase inhibitor, Sunitinib was reported to suppress various kinds of tumors, such as human metastatic renal cell carcinoma, lung carcinoma, glioblastoma, melanoma [12][13][14][15] by inhibiting the activation of VEGFR-1, 2 and 3, PDGFR-α, β, Kit, Flt3 and RET. It has been approved by FDA to treat advanced renal cell cancer and imatinib mesylate-intolerant or resistant gastrointestinal stromal tumor since 2006 [16,17].
Our results showed that R8 is more potent than Sunitinib in suppression of all types of cancer cells tested in this study, with much lower IC 50 values when compared to Sunitinib. Moreover, when specifically focused on lung cancer, one type of cancer with high morbidity and mortality, R8 exhibits stronger inhibition than Sunitinib on colony formation, cell cycle arrest and migration of lung cancer cells A549 and NCI-H226, although less effective than Sunitinib in inhibiting tube formation of endothelial cells. PI3K-Akt-mTOR has been identified as a key signal pathway in regulating various aspects such as cancer cell proliferation, migration, apoptosis and autophagy [18,19], and in this study we did observed that R8 induced significant down-regulation of p-AKT and total-Akt than Sunitinib did. Because of its key role in promoting cell cycle from G1 to S phase, down-regulated p-AKT resulted in cell cycle arrest at G0/G1 phase and proliferation inhibition.
Suppression of the Akt/mTOR pathway may also trigger autophagy and apoptosis [20][21][22]. Autophagy represents an important cell physiologic process in response to starvation, hypoxia, ER stress and oxidative stress. Autophagosomes engulf cellular organelles, such as mitochondria, ribosome, then fuse with lysosome to recycle energy used for cell survival during stress [23]. Generally, autophagy is considered as a way to protect cells from apoptosis [24][25][26]. Nonetheless, if stress over a safe threshold, autophagy can lead to cell death, called autophagic cell death or type II programmed cell death [27][28][29]. Indeed, several studies have confirmed the complex relationship between apoptosis and autophagy during the progress of cell treated with chemotherapeutic agents [30][31][32][33]. In this study we checked apoptosis and autophagy status of cancer cells treated with Sunitinib or R8 after 48 h. Because of down-regulation of PI3K-Akt-mTOR pathway, R8 induced significantly autophagy than Sunitinib in a dose-dependent manner marked by accumulation of LC3-II. Bcl-2 family proteins were activated to trigger apoptosis, resulted in up-regulation of Bax, cleavage of caspase-3 and PARP [34]. In the current study, we observed that suppression of R8-induced autophagy resulted in attenuated cleavage of caspase-3 and PARP, while the cleavage of caspase-3 and PARP were elevated on inhibition of Sunitinib-induced autophagy. The different sensitivities to 3-MA and CQ observed in this study were due to different mechanisms the two inhibitors involved in autophagy. Although both 3-MA and CQ are autophagy inhibitors, they block different steps of autophagy. 3-MA suppresses an early stage of autophagy by inhibiting the class III PI3K to block the production of PI3P, which is essential for early autophagosome formation [35]. CQ is a lysosomotropic agent which blocks the late stage of autophagy by elevating the lysosomal pH, inhibiting autophagosome fusion with lysosome and autophagosome degradation [36]. Taken together, these results suggested that Sunitinib-induced a pro-survival autophagy to evade cytotoxicity in lung cancer cells, similar to that in ovarian cancer [37], whereas R8-induced autophagy may have synergetic effect on cell death with apoptosis.
Our results also showed some differences between A549 and NCI-H226. NCI-H226, with more abundant receptor tyrosine kinases, showed more resistant to TKIs with higher IC 50 than A549, while A549, with less hyperactivated tyrosine kinase receptors, was more sensitive to R8. Interestingly, we found that the inhibition of EGFR phosphorylation in A549 was not only independent of the dose of R8, but even showed up-regulation of p-EGFR when treated with higher dose of R8. Whether this is just a feedback of self-protective response to the severe stress or by some other specific mechanisms need to be investigated in the future. Furthermore, the differences between the two cell lines on critical steps during autophagy may also be due to the different responses to RTK inhibition. For A549, EGFR is a prominent tyrosine kinase. Activated EGFR phosphorylates Beclin-1 at three different tyrosine residues, Y229, Y233 and Y352, provokes the association of Beclin-1 with Rubicon/Bcl-2, which results in reduced autophagy. Exposure of A549 to tyrosine kinase inhibitors restores the Beclin-1-VPS34 interaction and triggers an increase in autophagy, but it is independent of mTOR, modulation of the Beclin-1-Rubicon/Bcl-2 complex is dominant in the EGFR-mediated autophagy. In NCI-H226 cell line, PDGFR can not induce Beclin-1 phosphorylation, but can strongly activate autophagy by the PI3K-mTOR pathway [38][39][40]. Thus, the mechanism of autophagy regulation, by EGFR-Beclin-1 or PI3K-mTOR, may be an important factor to decide critical steps in cell fate.
Notably, we demonstrated in this work that R8 can significantly inhibit the phosphorylation of PDGFRβ and VEGFR2. Previous studies have shown the association between dysregulation of these two RTKs and chemoresistance: increased PDGFRβ and VEGFR2 protein levels were associated with resistance to platinumbased chemotherapy and adverse outcome of ovarian cancer patients [41]; increased VEGFR2 gene copy was associated with chemo-resistance and shorter survival in non-small cell lung cancer patients who received platinum adjuvant chemotherapy after surgical resection [42]. Moreover, pretreatment with VEGFR and PDGFRβ inhibitors can inhibit tumor growth in a mouse colorectal carcinomatosis model [43]. Therefore, malignancies with dysregulated PDGFRβ and VEGFR2 may benefit from R8 treatment due to the potent inhibition on these RTKs.
In conclusion, current work showed that the novel compound R8 is a potent multikinase inhibitor which can effectively inhibit the phosphorylation of PDGFRβ, VEGFR2, EGFR and C-Kit. R8 shows more in vitro suppression on lung cancer cells than Sunitinib does by the synergetic effect of autophagic cell death and apoptosis. This provides the potential that R8 alone may serve as a good chemotherapeutic drug without co-administration of an autophagy inhibition drug to enhance the therapeutic efficacy which is suggested for many chemotherapeutic drugs. Further in vivo assays such as non-small cell lung cancer patient-derived xenograft models which express high level of PDGFRβ, VEGFR2, EGFR, C-Kit will be carried out to demonstrate that whether R8 can be developed into a new member of MKI drug family.

Cell culture
Human gastric carcinoma cells SGC7901, non-small cell lung cancer cell line A549 and NCI-H226, laryngeal carcinoma cell line HEp-2, renal adenocarcinoma cell line 786-0, colorectal cancer cell line SW620 were obtained from American Type Culture Collection (Manassas, VA, USA). Cells were maintained in RPMI-1640 medium, supplemented with 10% fetal bovine serum (FBS) at 37°C in a humidified incubator with 5% CO 2 , Primary human umbilical vein endothelial cell line (HUVEC) was a kind gift from Professor Jianan Wang (Zhejiang University School of Medicine, China).

Cell viability assay
R8 and Sunitinib were dissolved in ddH 2 O or RPMI-1640 medium respectively, to give a storage concentration of 20 mM and stored at -20°C for no more than one week. For A549, cells were seeded at a density of 8×10 3 per well in 96-well plate and cultured for 24 h. After 24 h, 10% FBS-RPMI1640 was removed, cells were starved overnight. Then, cells were exposed to different concentration of R8 or Sunitinib (range from 2.5 μM-80 μM) for 1 h before stimulation. Subsequently, cells were stimulated with same volume of 1% FBS-RPMI1640 and cultured in the presence of R8 or Sunitinib. After 48 h, cell viability was assessed by cell counting kit-8 assay as described previously [44]. The IC 50 of R8 and Sunitinib were calculated based on the cell viability curve.

Cell cycle assay
A549 (4×10 5 /well) and NCI-H226 (5×10 5 /well) cells were seeded in 6-well plate and cultured for 24 h. All procedures were the same as cell viability assay except for that we harvested the cells at 24 h after stimulation. Then cells were fixed with 75% ethanol overnight followed by RNase treatment and propidium iodide (PI) staining as the instruction of PI/RNase Staining Buffer indicated. Cells were subjected to flow cytometer (FACSCalibur, BD, San Jose, CA, USA) and analyzed using CELL Quest 3.0 software (BD, NJ, USA).

Colony formation assay
A549 (500/well) and NCI-H226 (800/well) cells were seeded in 6-well plate and cultured for 24 h. After serum-starved overnight, cells were exposed to different concentrations of R8 and Sunitinib for 1 h then stimulated with same volume of 10% FBS-RPMI1640. After that, cells were cultured for 14 days to allow colony formation. Colonies were fixed with 4% PFA for 15 min then stained with 1% crystal violet for 15 min. Colony consisting of 100 or more cells was scored [45]. The results were photographed by ChemiDoc XRS+ image system of Bio-Rad (Hercules, CA, USA) and Lecia DMIL LED inverted microscope (Wetzlar, Germany). www.impactjournals.com/oncotarget Migration assay A549 (6×10 5 /well) and NCI-H226 (5×10 5 /well) were seeded in 6-well plate and cultured for 24 h. After that medium was replaced by RPMI-1640 medium without FBS, cells were starved for 24 h. Then on the third day, 10 μl pipette tip was used to scribe gaps in each well. Scratched cells were washed with PBS for three times. Simultaneously, the gaps were photographed, set as 0 h. Finally, different concentrations of R8 and Sunitinib were added into the medium. Photographs were captured at 12, 24, 36 h by inverted microscope for the same position [46]. The migrated distances were analyzed by Image-pro plus software.

Tube formation assay
Matrigel basement membrane matrix was placed in 4°C overnight before use. On the following day, a 96-well plate was coated with 50 μl Matrigel on ice and incubated at 37°C for 30 min. Then HUVECs (10 4 / well) were added in serum-free medium in the absence or presence of various concentrations of R8 or Sunitinib. Each conditional group contained 5 wells. Pictures of representative ×10 fields were taken from 4 to 6 h. The length of endothelial tubes was measured by Image-Pro Plus to assess the tube formation ability.

Western blot analysis
Cells were seeded in 100 mm culture dishes at a density of 1.5×10 6 per dish. To evaluate the function of R8 on inhibiting growth factor stimulated receptor phosphorylation, cells were exposed to R8 for 1 h then harvested immediately on ice after simulated with 10% FBS-RPMI1640 for 10 min. To evaluate autophagy and apoptosis, cells were exposed to R8 and Sunitinib for 1 h then stimulated with 1% FBS-RPMI1640. Cells were harvested after 48 h treatment. For autophagy inhibitor, cells were treated with 3-MA (10 mM) or CQ (20 μM) for 2 h before exposed to R8 and Sunitinib.
The cell pellets were lysed in RIPA lysis buffer containing 100× Phosphatase Inhibitor Cocktail 2 and 50× Complete ™ Protease Inhibitor Cocktail for 1 h on ice. Cell lysates were centrifuged at 14,000×g for 15 min at 4°C then supernatant extracts were collected. The protein concentration was determined by Dc Protein assay kit of Bio-Rad. Equal amount of proteins (30 μg/ well) were separated by electrophoresis on 10-12% SDS-PAGE and transferred to PVDF membranes from Bio-Rad (Hercules, CA, USA). Membranes were blocked with 10% BSA-TBST for 1 h at room temperature and then incubated with primary antibodies at 4°C overnight. Next, membranes were incubated with an HRP-conjugated secondary antibody for 1 h at room temperature. The immune-reactive bands were visualized by the ECL kit from Millipore (Billerica, MA, USA).

Statistical analysis
Each experiment was performed at least for three times, and the data were shown as means±standard deviation. The statistic difference was analyzed by twotailed Student's t-test. P values <0.05 were considered to be statistically significant.

Author contributions
YX, JC participated in research design; YX, CL conducted experiments; YH, ZJ contributed new reagents and analytic tools; YX performed data analysis; YX, JC contributed to the writing of the manuscript.


Set-Based Tests for Genetic Association Using the Generalized Berk-Jones Statistic

Studying the effects of groups of Single Nucleotide Polymorphisms (SNPs), as in a gene, genetic pathway, or network, can provide novel insight into complex diseases, above that which can be gleaned from studying SNPs individually. Common challenges in set-based genetic association testing include weak effect sizes, correlation between SNPs in a SNP-set, and scarcity of signals, with single-SNP effects often ranging from extremely sparse to moderately sparse in number. Motivated by these challenges, we propose the Generalized Berk-Jones (GBJ) test for the association between a SNP-set and outcome. The GBJ extends the Berk-Jones (BJ) statistic by accounting for correlation among SNPs, and it provides advantages over the Generalized Higher Criticism (GHC) test when signals in a SNP-set are moderately sparse. We also provide an analytic p-value calculation procedure for SNP-sets of any finite size. Using this p-value calculation, we illustrate that the rejection region for GBJ can be described as a compromise of those for BJ and GHC. We develop an omnibus statistic as well, and we show that this omnibus test is robust to the degree of signal sparsity. An additional advantage of our method is the ability to conduct inference using individual SNP summary statistics from a genome-wide association study. We evaluate the finite sample performance of the GBJ though simulation studies and application to gene-level association analysis of breast cancer risk.

Introduction
Genome-Wide Association Studies (GWAS) have been successful in identifying the associations between thousands of Single Nucleotide Polymorphisms (SNPs) and a variety of complex traits (Manolio et al. 2009). A traditional GWAS analysis tests for the effect of each individual SNP, and this approach has shown that single-SNP effects are often weak across the genome (Visscher et al. 2012). Recently, set-based tests which jointly analyze a group of SNPs -e.g. SNPs in a gene, pathway, or network -have become increasingly popular as complementary tools which can boost analysis power in GWAS (Wu et al. 2010). These tests are also standard for rare variant analysis in whole genome sequencing studies (Lee et al. 2014).
SNPs can be aggregated into sets based on a variety of genomic features. For example, they can be grouped by physical position, such as location in a gene or Linkage Disequilibrium (LD) block, or similar biological functions, such as membership in a genetic pathway or protein network. Set-based analyses then allow for some natural advantages over individual SNP methods. Besides reducing the number of multiple comparisons across the genome, SNP-set methods can increase power by pooling sparse and weak effects into a stronger aggregated signal, as well as by incorporating biological information into the test (Wu et al. 2011). In addition, set-based interpretations of association may be more meaningful than their single-marker counterparts, such as in a gene-level or pathway-level analysis.
A number of set-based tests for genetic association studies have been developed in recent years, including burden tests (Li and Leal 2008), the Sequence Kernel Association Test (SKAT) (Wu et al. 2011), the minimum p-value test (MinP) (Conneely and Boehnke 2007), and most recently the Generalized Higher Criticism (GHC) (Barnett et al. 2017). SKAT and burden tests are examples of methods more suitable for detecting dense signals. If the signals reside in only a few SNPs that are not correlated with noise SNPs, then the power of SKAT and burden tests will suffer.
While certain SNP-sets may contain a large number of signals, it is more common that genomic constructs formed with GWAS data will have only a few signal SNPs. Interestingly, within tests designed for this sparse alternative setting, there are still subtle differences in performance. Under extreme sparsity, as in the case of only one or very few signals in the entire set, the minimum p-value test and the GHC have good power for detecting a SNP-set effect. However GHC and MinP can lose power under moderate sparsity settings, which are relatively common in gene and pathway level analyses of GWAS data. For example, in the Cancer Genetic Markers of Susceptibility (CGEMS) GWAS for breast cancer risk (Hunter et al. 2007), four out of 42 SNPs in the FGFR2 gene, a known breast cancer risk loci, showed strong evidence of association without reaching genomewide significance. It is hence of substantial interest to develop testing procedures that can reliably detect associations across a range of alternatives in the sparse signal regime.
When factors in a set are independent, several goodness-of-fit type methods have been proposed to perform set-based tests in the presence of sparse signals (Donoho and Jin 2004;Jager and Wellner 2007;Walther 2013). These methods test for the effect of a set by aggregating evidence from marginal test statistics, and they have been shown to possess attractive asymptotic properties when the size of a set goes to infinity. Specifically, they reach a so-called detection boundary when signals are sparse. In a certain sense, they are able to detect the weakest signals detectable by any statistical procedure under the sparse alternative. The class of tests with this ability includes the Higher Criticism (HC) (Donoho and Jin 2004) and the Berk-Jones (BJ) (Berk and Jones 1979). Compared to the HC and BJ tests, the minimum p-value test, for example, is known to attain the detection boundary over a smaller portion of the sparse regime. In terms of finite sample performance, it has been demonstrated through simulation that Berk-Jones outperforms Higher Criticism over a range of moderately sparse alternatives when marginal test statistics are independent (Walther 2013;Li and Siegmund 2015). Donoho and Jin (2004) provide an explanation for this result by showing that HC disproportionately weights evidence from the most extreme observed marginal test statistic, at the cost of losing sensitivity to detect signals in other locations.
However, a direct application of BJ to SNP-set testing is not desirable, as standard p-value calculations for BJ rely on independence between elements in a set (Moscovich-Eiger and Nadler 2017). This assumption is violated by LD-induced correlation between neighboring SNPs in a SNP-set. In addition, we will see that even if we correct the inference procedure of the Berk-Jones, the power of BJ can fall dramatically in high LD settings.
To overcome the challenges posed by correlated SNPs, this paper proposes the Generalized Berk-Jones (GBJ) statistic for testing the association between a SNP-set and outcome. GBJ accounts for LD among SNPs in a set while still preserving the ability to detect moderately sparse and weak signals in finite samples. In fact, GBJ reduces to the Berk-Jones statistic when all SNPs in a set are independent. GBJ can also be applied to SNP-set tests using both individual-level genotype data or GWAS summary statistics from single SNP analysis. To facilitate use, we additionally provide an analytic p-value calculation for GBJ. Our method is more computationally efficient than permutation and is shown to be accurate even at the extremely small levels required for genome-wide significance.
Additional insight into the strengths and weaknesses of GBJ is provided by studying the rejection regions of SNP-set tests developed for sparse alternatives. The rejection regions allow us to quantitatively describe how the power of each test is susceptible to changes in parameters such as the amount of correlation between SNPs or the size of the SNP-set. Since in practice we never have knowledge of the type of alternative, we also propose an omnibus test that combines GBJ, GHC, MinP, and SKAT for added robustness to different degrees of sparsity. An extensive simulation study demonstrates that GBJ outperforms alternative methods in testing SNP-set effects when signals are weak and moderately sparse, and we also show that the omnibus test is robust to a wide range of sparsity levels.
The remainder of the paper is organized as follows. Section 2 discusses the SNPset testing framework using both individual-level data and GWAS summary statistics. In Section 3 we propose the Generalized Berk-Jones statistic for testing the association between a SNP-set and outcome. We also provide an analytic p-value calculation for GBJ and develop the omnibus test. Section 4 compares the rejection regions of GBJ and other tests designed for the sparse regime. In Section 5 we demonstrate the finite sample performance of GBJ through simulation. Section 6 presents an analysis of the CGEMS data, and we conclude with a discussion in Section 7.

Individual-Level Genotype Data
We begin by describing the SNP-set testing framework using individual-level data on genotype, outcome, and other covariates for n total subjects. Suppose for subject i, i = 1, ..., n, we observe the outcome Y i , a genotype vector G i = (G i1 , ..., G id ) T of d SNPs in a SNP-set, and a vector of q covariates X i = (X i1 , ..., X iq ) T . Assume that Y i conditional on (G i , X i ) follows a distribution in the exponential family (McCullagh and Nelder 1989) with the density function f , and c(·), are known functions, θ i is a canonical parameter, and φ is a dispersion parameter. Let is a differentiable monotone link function. We here only consider canonical link functions for simplicity. In matrix notation, the data take the form The null hypothesis of no association between a SNP-set and outcome, after controlling for covariates, is given by H 0 : β = 0. Both the size of the set d and the sparsity of signals can vary greatly between sets, e.g. from gene to gene, and the number of nonzero β j is unknown. Our aim is to develop a test suitable for different levels of sparsity while also accounting for the correlation among individual SNP test statistics.
The marginal score statistic for β j , j = 1, ..., d, is Under the null, the d marginal score test statistics have an asymptotic multivariate normal distribution where Σ jj = 1 for all j, and for j = k we can estimate (1)

GWAS Summary Statistics
Many GWAS may not release individual-level data due to logistical challenges or data confidentiality agreements. Instead it is much more likely that a marginal test statistic for association with the outcome is available for each individual SNP (Pasaniuc and Price 2016). It is hence of great interest to be able to perform SNP-set testing using precomputed Z j from across the genome. To test a set of precomputed Z j with GBJ, we require estimation of their correlation matrix Σ using external information.
Assume we have a panel of reference genotypes from n r subjects of the same ethnicity as those used to construct the summary statistics. For example, this could come from the publicly available 1000 Genomes dataset (1000Genomes Project Consortium 2015. We estimate Σ using equation (1) but replace G .j and X with G (r) .j and X (r) , where G (r) .j is the n r × 1 genotype vector of SNP j from the reference panel, X (r) = (1, PC 1 , ..., PC m ), PC 1 , ..., PC m are the first m principal component vectors calculated from the reference panel, and m is the same number of principal components as was used to control for population stratification (Price et al. 2006) in the original GWAS analysis of the data. Additionally for each subject we estimate v( µ 0i ) by setting µ 0i equal to the sample mean of the outcome. For a normally distributed outcome, this is exact as v(·) = 1. For a binary outcome, since population stratification is the primary confounder of the SNP-outcome relationship in GWAS, and because µ i0 generally varies slowly with the principal components, this approximation is practically reasonable. Ultimately we are approximating .k up to a scale parameter, with the scale parameter eventually cancelled out in Σ jk .

The Berk-Jones Statistic
We briefly review the Berk-Jones statistic in this section to help introduce the Generalized Berk-Jones statistic in Section 3.2. The BJ statistic is designed to test for H 0 : β = 0 against the alternative that a nonempty subset of the β j are nonzero, assuming the marginal test statistics are independent. LetΦ(t) = 1 − Φ(t) denote the survival function of a standard normal random variable and Φ −1 (t) denote its inverse. Let |Z| (j) denote the order statistics of the vector that results from applying the absolute value operator to each element of Z, so that |Z| (1) is the smallest value of Z in magnitude. Set which is the number of marginal test statistics with a magnitude greater than or equal to some threshold t. For a fixed t ≥ 0, and if Z j iid ∼ N (0, 1) for all j, then S(t) has a binomial distribution with size d and mean parameter π = 2Φ(t). This viewpoint motivates the Berk-Jones statistic for independent observations (Donoho and Jin 2004) as: and the second line uses the characterization of S(t) as a binomial random variable. We see that BJ can roughly be explained as the maximum of a one-sided likelihood ratio test on the mean parameter of S(t), performed over the larger half of observed test statistic magnitudes. At t = |Z| (d−j+1) we have π = 2Φ |Z| (d−j+1) under the binomial likelihood null, and we have π = j/d under the binomial likelihood alternative. We say binomial likelihood null and alternative to make clear that we are talking about an interpretation of the Berk-Jones statistic and to distinguish from the actual set-based null and alternative hypotheses being tested.
Note that the Higher Criticism test differs from the Berk-Jones by replacing the likelihood ratio statistic in (2) with the Pearson Chi-square statistic {S(t)−2dΦ(t)} 2 /{2dΦ(t)(1− 2Φ(t)}. Let k be the number of causal SNPs in a set. The sparse regime is designated as k < d 1/2 , and we call d 1/4 < k < d 1/2 moderately sparse, with k ≤ d 1/4 referred to as extremely sparse. Donoho and Jin (2004) showed that, when the Z j are all mutually independent, both HC and BJ are able to reach the detection boundary over the entire sparse signal regime as d → ∞. Walther (2013) and Li and Siegmund (2015) showed that the BJ statistic generally has better power than HC when the size of the set d is finite and the signals are moderately sparse.
If Z 1 , ..., Z d are correlated, as they will be for test statistics arising from neighboring SNPs in a gene, then S(t) no longer has a binomial distribution under the null. In this case, the standard Berk-Jones statistic no longer has a meaningful interpretation, and we may expect it to lose efficiency. In fact, we will show later that the rejection region of the Berk-Jones has a less desirable shape under various correlation structures, leading to a significant loss in power when the test statistics are not independent. Therefore we are interested in developing a modified BJ statistic that can account for correlation among the marginal test statistics in a set and thus possesses rejection regions which are more robust to arbitrary correlation structures.

The Generalized Berk-Jones Statistic
We now propose the Generalized Berk-Jones statistic for testing the association between a SNP-set and outcome. Following the spirit of Berk-Jones, GBJ considers a likelihood ratio type statistic on the mean parameter of S(t), but the key difference is GBJ explicitly accounts for the correlation structure of Z 1 , ..., Z d . More precisely, we define the GBJ statistic as: When the Z j are correlated, S(t) follows either an underdispersed or overdispersed binomial distribution instead of the standard binomial. However finding the exact distribution of S(t) when Σ = I is difficult. For a general Σ, computing Pr {S(t) = m} requires iterating through d choose m terms and is very time consuming. In special cases, such as when Σ has an exchangeable correlation structure with Σ = (1−ρ)I+ρ11 T , the calculation is much easier. However these scenarios occur rarely, if ever, in practice.
We propose to approximate the full distribution of S(t) using an Extended Beta-Binomial (EBB) distribution (Prentice 1986). The Extended Beta-Binomial is a reparameterization and extension of the standard Beta-Binomial(α, β) distribution with the standard Beta-Binomial being a special case of the EBB. A random variable V ∼ EBB (d, λ, γ) has probability mass function where we follow the convention a k=0 c k = 1 for a < 0. The mean of V is given by The Extended Beta-Binomial distribution reduces to the Beta-Binomial distribution if we set λ = α(α+β) −1 and γ = (α+β) −1 for α, β > 0. Because the standard Beta-Binomial distribution requires α, β > 0, it cannot accommodate underdispersion and never reduces to the binomial distribution. In contrast, the EBB allows for both overdispersion and underdispersion, and it reduces exactly to the binomial distribution when γ = 0. This mechanism allows our GBJ statistic to reduce to the Berk-Jones when there is no correlation among the observations.

Calculation of the Generalized Berk-Jones Statistic
We now describe more precisely the mechanics of calculating the GBJ statistic. To begin, check if the condition 1 2Φ |Z| (d−j+1) < j/d is satisfied at any j ≤ d/2. If the condition is never satisfied, then the observed value of GBJ is 0 and we do not need to perform any more computation. The following steps should only be taken on indices j ≤ d/2 where the condition is satisfied.
At each qualifying j, we approximate the distribution of S(|Z| (d−j+1) ) by an Extended Beta-Binomial random variable under both the binomial likelihood null and the binomial likelihood alternative. Denote these two variables by V , where E 0 and Var 0 denote the expectation and variance conditional on a ) using the same equations except with E a and Var a , which are the expectation and variance conditional on The first moment matching equation is simple to solve, since clearly The variance term in the second equation is more difficult. We can use Theorem 1 of Barnett et al. (2017) for Var 0 . For Var a , we need the following theorem: The proof of this theorem is given in the Supplementary Materials. The terms in the infinite sum shrink very quickly, and in practice we see good accuracy using only the first ten.
After matching all four parameters (λ The maximum value of GBJ (j) d among all qualifying j is then the observed Generalized Berk-Jones statistic.

Analytic P-value Calculation
Let G d be a general supremum-based global statistic such as the GBJ statistic. Suppose G d is constructed using independent marginal test statistics Z 1 , ..., Z d . Denote the observed value of this statistic by g, where higher values of g indicate more evidence for the alternative. As noted by Moscovich-Eiger and Nadler (2017), the p-value for g can often be written boundary points' that come from inversion of the test statistic. The points b 1 , b 2 , ..., b d will depend on g and d and will be different for different choices of a global statistic, but for the sake of presentation, we will suppress this dependency in the notation. Moscovich-Eiger and Nadler (2017) proposed a method that can calculate the p-value of G d very quickly if Z 1 , ..., Z d are independent. However when Z 1 , ..., Z d are correlated, their techniques for a fast calculation are not applicable.
An exact p-value for GBJ, and for any global test applied to correlated observations, must take into account the covariance structure of Z. The p-value for these tests is then where b j is understood to additionally depend on Σ as well. We are unaware of any computationally feasible expressions to calculate the joint distribution of the order statistics |Z| (1) , ..., |Z| (d) when d is moderate or large. The Supplementary Materials provides a procedure to compute this probability by partitioning the region into d! separate sections, but the method is very computationally intensive and not feasible for use with d > 10. However, an alternative way to write the rejection region of equation (5) is The right hand side of (6) suggests that the quantity can be calculated recursively. Indeed, The quantity in (6) is just q d,0 and can be calculated recursively as follows from steps similar to the proof of Theorem 1. Note that we can generalize the scheme described above to calculate p-values for many different supremum-based global tests by adopting the general approach of Moscovich-Eiger and Nadler (2017). As long as the test statistic can be inverted to create the bounds b 1 , ..., b d , we can use equation (7) to calculate its p-value when applied to correlated observations. In particular, we can use this procedure to perform p-value calculations for the HC, GHC, BJ, and GBJ statistics. Both the calculation and the Generalized Berk-Jones test are implemented in the R package GBJ, available on the CRAN repository.

The Omnibus Test
While we will show that the Generalized Berk-Jones test possesses an attractive finite sample rejection region when signals are moderately sparse, GBJ may also lose power in the presence of very sparse or dense signals. As SNP-set inference involves testing for a composite alternative H 1 : β = 0, there is no uniformly optimal test for both sparse and dense alternatives. Since signal sparsity varies between genes, the best test will also change from gene to gene, but it is unknown prior to scanning the genome. Thus we propose an omnibus test that offers robust power over a range of different sparsity levels.
The omnibus test is constructed by combining the SKAT, GBJ, GHC, and minimum p-value statistics, which have been described above. The motivation for choosing these four methods is to combine tests that are known to have good power when signals are dense, moderately sparse, very sparse, and the sparsest possible, respectively. The MinP method uses the set's largest marginal test statistic in magnitude |Z| (d) as a test statistic. When the Z j are independent, Donoho and Jin (2004) showed that MinP asymptotically reaches the same detection boundary as HC and BJ in the very sparse regime k ≤ d 1/4 but not the moderately sparse regime d 1/4 < k < d 1/2 . In finite samples, MinP can have better power than the other three tests when there are only one or two causal SNPs. In contrast, SKAT is known to have high power when signals in a SNP-set are dense.
The omnibus test first applies each of the four tests to the same SNP-set, and then it carries forward the smallest p-value from the four tests as a test statistic. Specifically, the omnibus test statistic is defined as: where p GBJ , p GHC , p SKAT , and p M inP denote the p-values of the four tests applied on the same SNP-set. As these tests are applied to the same data, the four p-values will be correlated.
Calculations of the p-value for OMNI must again account for the correlation between tests. We employ a Gaussian copula approximation for the joint distribution of the inverse-normal transformed p-values: where Φ M (·; R) denotes the joint cumulative distribution function of a multivariate normal distribution with mean vector zero and correlation matrix R. The correlation matrix R of the four component test statistics is estimated through parametric bootstrap under the null. For each subject i in the study, we simulate a new outcome based on the null mean µ 0i . When individual-level data are not available, we take µ 0i to be the same constant for all subjects as an approximation. Then each of the four tests are applied with the simulated outcome instead of the original one. The original design matrix, or approximated design matrix if working with summary statistics, is used each time. Each of the four p-values is subtracted from 1 and then inverse-normal transformed; under the null hypothesis the four transformed values have marginal normal distributions with mean zero. As we only need to estimate the correlation matrix R, only a small number of parametric bootstrap samples are needed. In practice, this procedure is repeated 100 times, and then we set R equal to the sample correlations of the inverse-normal transformed statistics. We will see that this omnibus test performs well across a variety of settings.

Rejection Region Analysis of Different SNP-Set Tests
We study in this section the finite sample rejection regions for the BJ, GBJ, HC, and GHC tests, and we advocate for viewing these statistics as boundary-defining algorithms. Consider a fixed SNP-set s with size d = d s and correlation structure Σ = Σ s , and suppose we wish to conduct inference at level α = 0.01. Using the p-value calculation from above, we can employ standard root-finding routines to find the observed value g s which would result in a GBJ p-value of 0.01 for a SNP-set with these characteristics. Then inverting g s to find boundary points b 1 , b 2 , ..., b ds as in Section 3.4 constructs a rejection region in terms of in terms of |Z| (1) , ..., |Z| (ds) . That is, if the observed value of |Z| (j) were larger than b j for any j = 1, 2, ..., d s , then the GBJ p-value for this set would be less than 0.01. Plotting the bounds for various tests using different values of α, d, and Σ shows us exactly what types of signals a given test is well-powered to detect at level α. For the same setting, a test with smaller bounds is preferred, as it will provide more finite sample power.
To numerically compare the rejection boundaries, consider a simplified model of SNP-set correlation structure where the set is partitioned into only two sections. Let one section be the independence section, where all SNPs in this portion are completely independent of all other SNPs in the set. Let the other section be the correlated section, where all SNPs in this portion have common pairwise correlation ρ with other SNPs in the section. For our numerical study, ρ = 0.3 for the correlated section. We investigate SNP-sets of size d = 20 and 100, correlated sections which contain 50% and 75% of the SNPs, and tests at size α = 0.01. These parameters are chosen to represent reasonable boundaries on the correlation structures seen in common GWAS data; Dawson et al. (2002) estimated that the average r 2 between SNPs separated by 100kb is around 0.1.
The rejection regions for each SNP-set are plotted in Figure 1. At the jth coordinate on the x-axis, if the observed |Z| (j) lies above the boundary of a particular test at that coordinate, then we would reject the null for that test at level α = 0.01. The lines on the graph are added to aid in visualization, but there should be no interpretation of interpolation between the points. It does not make sense to think of the boundary at |Z| (2.5) , for example. While standard methods for inference on HC and BJ are invalid in the presence of correlation, valid p-values for these tests can be computed with the same ideas we have introduced for GBJ inference, specifically following equations (5)-(7). Thus we can show that HC and BJ sometimes have much less desirable rejection regions when SNPs in a set are correlated.
One of the clearest trends from Figure 1 is that the HC and GHC boundaries are lower for a small region around |Z| (d) , and then the BJ and GBJ boundaries quickly become smaller as we move left. This behavior indicates that HC and GHC are better at detecting the sparsest alternatives with only a very few signals, as those signals would almost always manifest as the test statistics with the largest magnitude. In contrast, the plots demonstrate that BJ and GBJ can have more power to detect weaker, less sparse signals which may be more easily found by examining the test statistic which is, say, fifth or tenth largest in magnitude. The boundaries of HC and GHC can drop below BJ and GBJ again for the smallest observed magnitudes, but signals would only be found in these observations if they are particularly dense, a setting which is not the focus of our efforts. The intuition we can glean from this figure is closely aligned with the theoretical development of Donoho and Jin (2004) and the simulations of Li and Siegmund (2015) when the marginal test statistics Z j are independent. These authors showed that HC is attuned to detect sparse signals arising at the very tail of the observed distribution, while BJ has more power as the number of signals rises.
These results also show why BJ is likely to have low power for detecting sparse signals when the level of correlation is high. When 75% of the SNPs are correlated, the rejection boundary for BJ at the largest few observations is the highest by multiple orders of magnitude on the p-value scale. It would not be desirable to apply BJ in these types of settings, as the test loses an extremely large amount of sensitivity to detect signals in the most outlying values. BJ is still likely to be suitable for detecting dense signals in these situations. Here, GBJ acts as a compromise between BJ and GHC under high correlation.  If the jth smallest test statistic in magnitude is greater than the boundary point for a given test at any point j on the x-axis, then we would reject the null using that test at level α = 0.01. The difference between BJ and GBJ becomes much more pronounced as both the size of the set and the amount of correlation increase.
GBJ provides a much lower boundary than BJ at the tail in exchange for slightly higher boundaries near the middle. Thus, GBJ can detect both sparse and dense signals in this example. On the other hand, GBJ provides a slightly higher boundary than GHC at the tail in exchange for lower boundaries past the tail, so it trades some power in the extremely sparse regime for more power to detect moderately sparse signals. We see that choosing a different statistic is essentially choosing a different boundarysetting algorithm, and this choice should ideally be informed by parameters such as the amount of correlation and estimated sparsity level. Ultimately these plots illustrate that there is no single best global test for all types of alternatives. A genome-wide analysis strategy using the omnibus test will be likely to have robust power across different sparsity settings, correlation structures, and SNP-set sizes.

Type I Error of the Generalized Berk-Jones Test
We first illustrate that our p-value calculations for the GBJ and omnibus tests are accurate enough to control the Type I error rate at levels required to declare genomewide significance of a SNP-set. To replicate the setting of traditional GWAS data, we perform the size simulation on random regions across chromosome 5 which correspond to known gene sizes. We also conduct the size simulation on a high-LD subset of the FGFR2 gene and a low-LD subset of the FGFR2 gene to parse LD-related effects. We choose FGFR2 because it contains both high and low LD regions, and because it will later be the most significant gene in our analysis of the CGEMS data. All SNP-sets contain genotypes simulated with HAPGEN2 (Su et al. 2011) using the CEU population from HapMap3 as a reference panel.
In all simulations the outcome is generated as Y ∼ N (0, 1), and we fit the linear regression model (1) with β = 0 and X i = 1. Each simulation is repeated 20 million times, and we report the Type I error down to 10 −5 . Table 1 shows that our GBJ p-value calculation is accurate and protects the correct size for correlation structures seen in actual data. The p-value calculation for the omnibus test is similarly accurate at the most stringent significance levels, but it is somewhat conservative at larger significance levels. Table 1: Type I error of GBJ computed over 20 million simulations. The strong LD setting denoted by SLD refers to eight SNPs from FGFR2 which are highly correlated. The weak LD setting denoted by WLD refers to eight SNPs from FGFR2 which demonstrate only a small amount of correlation with each other. The Chr5 setting refers to over 500 regions on chromosome 5 corresponding to known gene sizes.

Power of GBJ Under Varying Hypothetical Correlation Structures and Sparsity Settings
To study the power of GBJ, we first conduct simulations under a variety of hypothetical correlation structures and sparsity settings. The performance of the GBJ is compared to the minimum p-value test, GHC, SKAT, and the omnibus test described in Section 3.5. The MinP test p-value is calculated by casting MinP as a boundary-defining test with b j = |Z| (d) for all j, and then computation proceeds through the methods described in Section 3.4. GBJ and GHC similarly use the p-value calculation described above. For SKAT, we run the corresponding R package.
To study how power is impacted by different correlation structures between the SNPs, we utilize block correlation structures which are slightly more complex than those used for the rejection region analysis in Section 4. Specifically, consider a set of causal SNPs that are correlated amongst themselves with common pairwise correlation ρ 1 . All other SNPs are then non-causal, and we allow half of them to have an exchangeable correlation structure with correlation ρ 3 ; the other half of the non-causal SNPs are completely independent of all other non-causal SNPs. Finally the pairwise correlation between a causal SNP and a non-causal SNP is set at ρ 2 . The three correlations ρ 1 , ρ 2 , ρ 3 will vary between 0 and 0.3. All SNPs are generated to have minor allele frequency of 0.3.
We demonstrate the effects of signal sparsity by using a large SNP-set of d=100 SNPs and varying the number of causal SNPs from k = 1 to k = 10. This allows us to examine power profiles in the very sparse regime (one to three causal SNPs), in the moderately sparse regime (four to nine causal SNPs), and at the edge of the dense regime (ten or more causal SNPs). The true disease model is where all the β j are the same and depend on the number of causal SNPs k. We reduce β j slightly as k increases in order to keep the power of each test below one throughout the entire sparse regime. The full details on effect size for each simulation are available in the Supplementary Materials. Figure 2 considers the case where the noise SNPs are independent, and Figure 3 considers the case where the noise SNPs are correlated. We perform 500 simulations at each different value of k and test at α = 0.01. All the power curves are smoothed to show empirical power. The first significant trend appearing in Figure 2 is the effect of sparsity on power. Among the non-omnibus tests, we see that GHC and MinP perform well when the number of causal SNPs is low, as these tests often have the most power in the very sparse regime. In both panels of Figure 2, the transition to GBJ having more power than GHC and minP occurs in the moderately sparse regime. Then as the number of causal SNPs increases into the dense regime, SKAT begins to catch up and eventually becomes the most powerful test. This behavior matches our intuition as well as previously published simulation results. GHC and minimum p-value place excess weight on the most outlying observations, so they are well-tuned to detect the very sparse signals. The rejection region of GBJ is better-suited to find moderately sparse signals, and SKAT is known to perform well with dense signals.
The relationship between sparsity and power can be modified by the total amount of correlation. In the left panel of Figure 3 we set ρ 1 = ρ 3 = 0.3, which corresponds to the situation where causal SNPs and non-causal SNPs are correlated within themselves, but the two groups are independent of each other. In this case, MinP and GHC become the top-performing tests for a larger range of sparsity settings, with GBJ losing some of its advantage in the moderately sparse regime. SKAT has almost no power in these situations, as the signals are sparse and there is no correlation between causal and noncausal SNPs. It appears that a large amount of correlation between the non-causal SNPs is detrimental to the performance of GBJ. An explanation for this behavior can be found . The correlation structure on the right is slightly simpler than the previous three structures, as we switch to an exchangeable correlation matrix in order to accommodate ρ 2 = 0.3 while still ensuring the correlation matrix is positive definite. As the number of causal SNPs increases, we decrease the effect size to keep power bounded away from one.
in the rejection region analysis of Figure 1. We see that the bounds of GBJ appear less favorable compared to GHC when the amount of correlation is high. Since over half of the SNPs in Figure 3 are correlated with ρ 1 = ρ 3 = 0.3, these settings represent a much larger amount of total correlation than was present in Figure 2. In the right panel of Figure 3 we investigate the setting of ρ 2 = 0 by using an exchangeable correlation structure, and SKAT dominates as the most powerful test across almost all sparsity levels. Here we break slightly from the above framework by using exchangeable correlation to accommodate ρ 2 = 0.3 while still ensuring the correlation matrix of the SNPs is positive definite. GBJ is a close second to SKAT under most sparsity settings with these parameters. SKAT is known to have good performance in the presence of LD between causal SNPs and noise SNPs, which makes signals appear to be dense. The increased density of signals also buoys the performance of GBJ compared to GHC and minimum p-value, which perform the worst under exchangeable correlation.
Never losing too much power compared to the best test, the omnibus test appears to be robust to LD structure and sparsity. This behavior is expected as our omnibus approach integrates information from tests which perform well across multiple sparsity settings. Thus we would anticipate that OMNI is more resilient than any single test.
GBJ also demonstrates good power in a variety of situations, and its overall strength and robustness across the entire sparse regime are attractive properties. In particular, GBJ is the best-performing test when the level of sparsity is moderate and there is weak correlation among the noise SNPs. GBJ is disadvantaged against minP and GHC when signals are extremely sparse or there is excess correlation among noise SNPs, but it outperforms these tests as signals become more dense. In contrast, GBJ provides slightly less power than SKAT as signals reach the dense regime or when there is moderate correlation between causal and non-causal SNPs, but GBJ is also much more robust than SKAT when signals are sparse and not correlated with noise SNPs. Similar to OMNI, GBJ does not fall far behind the best test in any given situation, which suggests that GBJ is a good choice to use when the signal sparsity is unknown. Power for the standard BJ is not plotted in the interest of space, but it behaves like a dense test, similar to SKAT, under correlation. This behavior again matches Figure 1, which showed that BJ is more suited to detect dense signals as the amount of correlation increases.

Power of GBJ under Actual Chromosome 5 Correlation Structures
We conduct one final simulation to investigate the power of Generalized Berk-Jones under the unstructured LD patterns found in real GWAS data. In this simulation, we choose blocks of 40 SNPs at random locations on chromosome 5, and then genotype data are generated using HAPGEN2. We choose 40 to again approximately match the characteristics of FGFR2. There are 2000 blocks chosen for each sparsity level, and 10 simulations are performed on each block, for a total of 20000 at each number of causal SNPs. Again, the effect size is decreased as the number of causal SNPs increases, with the outcome still generated from equation (8). Testing is performed at α = 1 · 10 −5 to mimic a practical analysis.
We see in the left panel of Figure 4 that GBJ, GHC, and the omnibus test all have very similar power curves in this setting, while SKAT and minP lag slightly behind. As the number of causal SNPs increases, GBJ demonstrates the best power by a small amount. These results are rather homogenous because sparsity levels are more coarse Recall that GBJ extended its advantage over GHC when there was less correlation among noise SNPs. When restricting our analysis to the blocks which have median |ρ 3 | < 0.1, we see a more substantial superiority of GBJ over GHC in the moderately sparse regime, similar to Figure 2. Thus it is possible to recover the patterns of the structured simulation in real genotype data. Obviously median |ρ 3 | is not a perfect summary measure, as it cannot single-handedly capture all the parameters in a 40 × 40 correlation matrix. Further parsing of the data would be necessary to see larger differences in performance. In a practical setting, we might switch between tests based on certain SNP-set characteristics, such as applying GBJ when the set is large and likely to have moderately sparse signals. These results do again demonstrate the robustness of GBJ across multiple situations, as it provides the most power across a large portion of the sparse regime.

Gene-Level Analysis of the CGEMS GWAS Data
The CGEMS breast cancer dataset contains a case-control sample of 1145 breast cancer cases, all postmenopausal women with European ancestry, and 1142 controls recruited from the Nurses' Health Study. These women were genotyped at approximately 550000 SNPs with the Illumina HumanHap500 array. The dataset was originally analyzed by Hunter et al. (2007) in the single-marker GWAS approach. The authors did not find any individual SNPs to reach the genome-wide significance level of 5 × 10 −8 , but they highlighted FGFR2 as a strong candidate for future studies based on four SNPs in the gene that showed suggestive association with breast cancer. Such a situation succinctly illustrates the burden of adjusting for multiple comparisons when testing individual SNPs. Gene-level analysis provides an attractive alternative strategy that can reduce the num-ber of comparisons and also aggregate evidence of signals across multiple SNPs in a gene. Here we perform a gene-level analysis to study the association between genes and breast cancer risk.
Since individual-level genotype data were available for this study, we first calculated the marginal test statistics for each SNP using the model in Section 2.1. Specifically, we fit a logistic regression model with four covariates -age and the first three genotype principal components to control for population structure (Price et al. 2006). Then, for each of 14991 genes, we collected the marginal test statistics for all SNPs located within the region defined by that gene. Each gene with more than one marginal SNP test statistic was analyzed with GBJ, GHC, SKAT, MinP, and the omnibus test.
In Table 2, we rank the top ten genes according to the smallest p-value produced by any of the five tests. In this sample, GBJ provides the strongest evidence of association for the top four genes and five of the top ten. Most of these genes are ranked highly by multiple other methods, however no other method provides the lowest p-value for more than two of the top ten genes. In fact, GHC and MinP produce the smallest p-value only once between the two of them. One possible explanation for the underperformance of GHC and MinP is that there may be multiple tagged SNPs surrounding the true causal loci for each of these genes, which could create a lack of extremely sparse alternatives. Table 2: Top significant genes in gene-level analysis of CGEMS breast cancer GWAS data, ranked by minimum p-value produced by any of the five tests. The test which produces the smallest p-value for each gene is highlighted in red.
Besides FGFR2, genes such as PTCD3 and POLR1A have also been implicated as risk loci in independent investigations (Boehm et al. 2007;Jia et al. 2011). The overlap of our findings with other studies and other statistics provides a level of reassurance that GBJ performs well in identifying truly significant genes and not simply spurious associations. Alternately, ABCA1 is an example of a gene that may not have received further scrutiny if we were not utilizing the GBJ test. ABCA1 expression has been linked with breast cancer risk (Smith and Land 2012), but MinP and GHC do not provide the same strength of evidence that GBJ does. It seems likely that there are more than a few signal SNPs in ABCA1, especially since ABCA1 contains a relatively large number of SNPs compared to the other genes in this dataset.
Perhaps due to the limited sample size, no test produces a p-value low enough to be declared significant after Bonferroni correction for 14991 genes. Still, this analysis highlights the advantages of Generalized Berk-Jones compared to alternative tests. The GBJ p-value for FGFR2 does come very close to the Bonferroni-corrected level (3.34 × 10 −6 ), and it certainly provides more evidence of association than the single SNP statistics. Additionally, GBJ often gives the highest measure of significance, and never the lowest, in the genes displayed, demonstrating its robustness across different set sizes and LD patterns.

Discussion
We have proposed the Generalized Berk-Jones statistic to test for association between a SNP-set and an outcome. Our GBJ generalizes the standard Berk-Jones by modifying the BJ statistic to directly account for the correlation between individual SNPs. This modification results in a test that is more powerful when SNPs are in LD. We also provide an analytic p-value calculation for GBJ and generalize it to a class of supremumbased global tests, allowing valid inference for HC, GHC, BJ, and other methods when these procedures are applied as SNP-set tests using correlated marginal test statistics. Rejection region analysis demonstrates that GBJ can be described as a compromise between Berk-Jones and Higher Criticism-type tests in terms of finite sample performance.
While our numerical analysis shows situations where GBJ does not set the lowest boundary at either |Z| (d) or |Z| (d/2) , GBJ generally comes very close to the lowest boundary at both locations, which affords it both robustness to signal sparsity and power to detect moderately sparse signals. GHC and HC often set the lowest boundary around |Z| (d) , but in return they concede a large amount of volume past the first few most extreme observations, which lowers power in the moderately sparse regime. BJ frequently sets the lowest boundary past the tail, but its tail boundary can be orders of magnitude larger than that of GBJ, HC, and GHC. Bounds in the expected signal regions must be viewed holistically, so slightly lower bounds at a few locations are not necessarily desirable if the price is much higher bounds in other signal locations, as in the case of BJ. Thus GBJ offers good power to detect moderately sparse effects without losing too much power when single-SNP signals are extremely sparse.
Simulation results reinforce the conclusions we find from examining the rejection regions of GHC and GBJ. Additionally we see that the MinP test performance is quite good when signals are very sparse, similar to GHC, but MinP does not perform as well as GHC when signals become more dense. SKAT has a unique power profile, as it can be particularly powerful when signals are dense or there is correlation between causal and non-causal SNPs, but it is also not robust to different correlation structures and will often have very little power in sparse settings where there is no correlation between causal and non-causal SNPs. The omnibus test offers robust power across different sparsity levels, and while it is rarely the best test, it also never has the worst power. When applied to data from the CGEMS study, we see that GBJ often produces the most significant p-values, perhaps owing to its versatility across different parameter settings.
In demonstrating that the BJ statistic can be adapted for increased robustness to correlation, we have also demonstrated that these types of boundary-defining algorithms can be modified to increase finite sample power under specific set-level parameters. It would be of interest to develop different boundary-defining methods that offer more favorable rejection regions in narrow but well-defined settings. For example, it may be possible to define tests which outperform GHC or minP over finer partitions of the extremely sparse regime.
In a similar vein, it would be interesting to understand the boundary shapes for other previously proposed boundary algorithms (Jager and Wellner 2007) in the class of Berk-Jones and Higher Criticism. While many of these algorithms share the same asymptotic guarantees of BJ and HC, little is known about their comparative finite sample performance, especially when observations in a set are correlated. These other methods might also have great value in specific settings as mentioned above. It would additionally be very valuable to develop SNP-set tests that are optimal in certain senses for arbitrary sparsity and correlation.
As genomic data collection techniques continue to evolve, it may be necessary to adapt the GBJ as well. In particular, the rise of whole genome sequencing and fine mapping studies is leading to the discovery of more SNPs with extremely rare minor alleles. Marginal test statistics generated from these SNPs are known to be non-Gaussian in finite samples, and thus they will not have the distribution we assume for GBJ. GBJ will need to account for null distributions that are not standard normal before SNP-sets containing rare variants can be tested.

SUPPLEMENTARY MATERIAL
The supplementary materials provide the proof of Theorem 1 from Section 3.3, offer further details on how to calculate the exact p-value from equation (5) in Section 3.4, and list the exact simulation parameters from Section 5.